Matlab 曲線擬合之 polyfit 、polyval、poly2str 函式
原文出處（僅供參考）
ttps://www.cnblogs.com/farewell-farewell/p/7227516.html
https://www.cnblogs.com/clairvoyant/p/4710015.html
1 Matlab 曲線擬合之polyfit與polyval函式
p=polyfit(x,y,n)
[p,s]= polyfit(x,y,n)
說明：x,y為資料點，n為多項式階數，返回p為冪次從高到低的多項式係數向量p。x必須是單調的。矩陣s用於生成預測值的誤差估計。
多項式曲線求值函式：polyval( )
呼叫格式： y=polyval(p,x)
[y,DELTA]=polyval(p,x,s)
說明：y=polyval(p,x)為返回對應自變數x在給定係數P的多項式的值。
[y,DELTA]=polyval(p,x,s) 使用polyfit函式的選項輸出s得出誤差估計Y DELTA。它假設polyfit函式資料輸入的誤差是獨立正態的，並且方差為常數。則Y DELTA將至少包含50%的預測值。
有如下資料
| 時間t |
1900 |
1910 |
1920 |
1930 |
1940 |
1950 |
1960 |
1970 |
1980 |
1990 |
2000 |
| 人口y |
76 |
92 |
106 |
123 |
132 |
151 |
179 |
203 |
227 |
250 |
281 |
1. y與t的經驗公式為 y = at^2 + bt + c
clear;
clf; %清除當前視窗
clc;
t = 1900:10:2000; %時間t
y = [76 92 106 123 132 151 179 203 227 250 281]; %人口y
plot(t,y,'k*');
hold on;
% figure; %重新開一個圖
p1 = polyfit(t,y,2);
plot(t, polyval(p1, t));
axis([1900 2000 0 300]); %影象xy軸範圍
disp(char(['y=',poly2str(p1,'t')],['a=',num2str(p1(1)),' b=',...
num2str(p1(2)),' c=',num2str(p1(3))]));
結果如下：
y= 0.0094289 t^2 - 34.7482 t + 32061.5711 a=0.0094289 b=-34.7482 c=32061.5711
2. y與t的經驗公式為y = a e^(bt)
clear;
clf; %清除當前視窗
clc;
t = 1900:10:2000; %時間t
y = [76 92 106 123 132 151 179 203 227 250 281]; %人口y
yy = log(y); %指數基尼必需的線性化變形
p2 = polyfit(t,yy,1);
b = p2(1);
a = exp(p2(2));
y2 = a * exp(b*t); %指數擬合函式式
plot(t,y,'rp',t,y2,'k-');
grid off;
xlabel('時間t');
ylabel('人口數(百萬)');
title('人口資料');
2 str2poly()和poly2str()
多項式表示法
MATLAB中採用行向量來表示多項式，將多項式的係數按降次冪次序存放於行向量中。即多項式P(x)=a0xn+a1xn-1+ ... +an-1x+an的係數行向量為P=[a0 a1 ... an],順序為高次冪到低次冪。舉個例子，x5-7x3+2x+4，可以用係數行向量[1 0 -7 0 2 4]來表示，多項式中缺少的次冪要用“0”補齊。
1、str2poly():
下面為str2poly()函式，實現從多項式的字串表示轉換為多項式的行向量表示，輸入不必降冪排列，程式碼如下：
1 %sr2poly.m
2 %把多項式的字串表示轉換為行向量
3 function Y=str2poly(X)
4 %X是字串形式的多項式
5 %Y是行向量形式的多項式
6 %輸入格式檢查
7 if(ischar(X)==0),
8 disp('輸入錯誤：輸入X必須是一個字串');
9 end;
10 %用正則表示式尋找'+'或者'-'的下標位置
11 index=regexp(X,'\+|\-');
12 %多形式的項數
13 L=length(index);
14 %用於儲存多項式的每一項資訊的單元字串矩陣
15 term=cell(1,L+1);
16 term(1)=cellstr(X(1:(index(1)-1)));
17 for i=1:L-1,
18 term(i+1)=cellstr(X(index(i):(index(i+1)-1)));
19 end;
20 term(L+1)=cellstr(X(index(L):end));
21 %如果第一項為空則刪除第一項
22 if(isempty(char(term(1)))),
23 term(1)=[];
24 %多項式的項數減一
25 L=L-1;
26 end;
27 %多項式係數矩陣
28 coefficient=[];
29 %多項式冪次矩陣，它與多項式係數矩陣一一對應
30 power=[];
31 for i=1:L+1,
32 %單項多項式字串表示式
33 substring=char(term(i));
34 %用正則表示式，尋找字串'x^',由於'^'是正則表示式中特殊字元，多以用'\^'
35 index2=regexp(substring,'x\^');
36 if(isempty(index2)==0),
37 %如果匹配上
38 if((index2(1)==1))
39 %單項多項式字串為'x^*'形式
40 coefficient=[coefficient 1];
41 power=[power str2num(substring((index2(1)+2):end))];
42 elseif(index2(1)==2)
43 %單項多項式字串為'+x^*'或者'-x^*','3x^*'形式
44 if(substring(1)=='+'),
45 coefficient=[coefficient 1];
46 power=[power str2num(substring((index2(1)+2):end))];
47 elseif(substring(1)=='-'),
48 cofficient=[coefficient -1];
49 power=[power str2num(substring(index2(1)+2:end))];
50 else 51 coefficient=[coefficient str2num(substring(1:(index2(1)-1)))];
52 power=[power str2num(substring((index2+2):end))];
53 end;
54 else
55 %單項多項式字串為'2.2x^'
56 coefficient=[coefficient str2num(substring(1:(index2(1)-1)))];
57 power=[power str2num(substring((index2+2):end))];
58 end;
59 else
60 %單項多項式字串不含'x^'
61 %用正則表示式，尋找字串'x'
62 index2=regexp(substring,'x');
63 if(isempty(index2)==0),
64 %如果匹配上'x'
65 if(index2(1)==1),
66 %單項多項式字串為'x*'形式
67 coefficient=[coefficient 1];
68 power=[power 1];
69 elseif(index2(1)==2),
70 %單項多項式字串為'+x'或者'-x','3x*'形式
71 if((substring(1)=='+')==1),
72 coefficient=[coefficient 1];
73 power=[power 1];
74 elseif(substring(1)=='-'),
75 coefficient=[coefficient -1];
76 power=[power 1];
77 else
78 %單項多項式字串為'2.2x*'
79 coefficient=[coefficient str2num(substring(1:(index2-1)))];
80 power=[power 1];
81 end;
82 else
83 coefficient=[coefficient str2num(substring(1:(index2-1)))];
84 power=[power 1];
85 end;
86 else
87 %單項多項式字串不含'x^','x',則斷定他是常數項
88 coefficient=[coefficient str2num(substring)];
89 power=[power 0];
90 end;
91 end;
92 end;
93 %合併同類項
94 N=max(power)+1;
95 Y=zeros(1,N);
96 for i=1:N,
97 index3=power==(N-i);
98 Y(i)=sum(coefficient(index3));
99 end;
在命令列輸入以下程式碼：(輸入時為字串形式)
1 >> str2poly('x^7-2x^6+3x^5-5x^4+x^2-1')
輸出程式碼如下：
1 ans = 2 3 1 -2 3 -5 0 1 0 -1
2、poly2str():
接下來是poly2str()函式，實現了把一個行向量表示的多項式轉換為常見的字串表示的多項式，程式碼如下：
1 %poly2str.m
2 %把多項式的行向量表示轉換為字串表示
3 function Y=poly2str(X)
4 %X是表示一個多項式的向量
5 %Y多項式的字串表示
6 %輸入檢查，如果X不是一個向量則退出
7 if isvector(X)==0,
8 disp('輸入錯誤：輸入X不是一個向量，請輸入一個代表多項式的向量！');
9 return;
10 end;
11 Y=''; %輸入字串
12 n=length(X);
13 for i=1:n, %把多項式的每一次冪轉換成字串
14 if(i~=1&&X(i)>0)
15 Y=[Y '+']; %若為正係數，必須新增‘+'
16 end;
17 %輸出係數
18 if(X(i)==0), %如果該次冪係數為零，則不輸出字串
19 continue;
20 elseif(X(i)==1&&i~=n), %如果該次冪係數為1，可以不輸出係數，只輸出想x^n
21 Y=Y;
22 else
23 Y=[Y num2str(X(i))]; %其他輸入情況
24 end;
25
26 if(i==n-1), %輸入x^n 1次冪輸入字串'x'
27 Y=[Y 'x'];
28 elseif(i==n), %0次冪不輸出x^n
29 Y=Y;
30 else
31 Y=[Y 'x^' num2str(n-i)]; %其他情況輸出x^n
32 end;
33 %如果不是最後一項，輸出'+'
34 end;
35 if(Y(1)=='+'), %修正如果0次冪為0時，造成末尾有多餘的字串‘+’
36 Y(1)=[];
37 end;
在命令列輸入以下程式碼：(輸入時為矩陣行向量形式)
1 >> poly2str([1 -2 3 -5 0 1 0 -1])
輸出程式碼如下：
ans =
x^7-2x^6+3x^5-5x^4+x^2-1