實驗目的：
1.Matlab中多項式的表示及多項式運算
2.用Matlab實現拉格朗日及牛頓插值法
3.用多項式插值法擬合數據
實驗要求：
1.掌握多項式的表示和運算
2.拉格朗日插值法的實現（參見呂同富版教材）
3.牛頓插值法的實現（參見呂同富版教材）
實驗內容：
1.多項式的表達式和創建；多項式的四則運算、導數與積分。
2.用Matlab實現拉格朗日及牛頓插值法。
3.用多項式插值法擬合數據。
實驗步驟：
1.多項式的表達式，MATLAB中使用以為向量來表示多項式，將多項式的系數按照降冪次序存放在向量中。多項式P(x)的具體表示方法：
的系數構成向量為：
。示例如下：
將向量表示的多項式用字符串輸出的通用函數示例：
例子
運行示例：
多項式的加法：
結果是
多項式乘法：
結果是
多項式除法：
多項式導數：
2.用Matlab實現拉格朗日，拉格朗日代碼：
1 function yi=Lagrange(x,y,xi) 2 m=length(x);n=length(y);p=length(xi); 3 if m~=n 4 error('向量x與y的長度必須一致'); 5 end 6 s=0; 7 for k=1:n 8 t=ones(1,p); 9 for j=1:n 10 if j~=k 11 t=t.*(xi-x(j))./(x(k)-x(j)); 12 end 13 end 14 s=s+t.*y(k); 15 end 16 yi=s; 17 end
運行示例：
牛頓插值法代碼：
1 function yi=newtonint(x,y,xi) 2 m=length(x);n=length(y); 3 if m~=n 4 error('向量x與y的長度必須一致'); 5 end 6 A=zeros(n); 7 A(:,1)=y; 8 for j=2:n%j為列標 9 for i=1:(n-j+1) %i為行標 10 A(i,j)=(A(i+1,j-1)-A(i,j-1))/(x(i+j-1)-x(i));%計算差商表 11 end 12 end 13 %根據差商表,求對應的牛頓插值多項式在x=xi處的值yi 14 N(1)=A(1,1); 15 for j=2:n 16 T=1; 17 for i=1:j-1 18 T=T*(xi-x(i)); 19 end 20 N(j)=A(1,j)*T; 21 end 22 yi=sum(N); %將x=xi帶入牛頓插值多項式，得到的yi的值 23 %A 輸出差商表 24 end
運行實例：
等距節點的牛頓向后插值代碼：
1 function yi=newtonint1(x,y,xi) 2 h=x(2)-x(1);t=(xi-x(1))/h; 3 n=length(y);Y=zeros(n);Y(:,1)=y'; 4 for k=1:n-1 5 Y(:,k+1)=[diff(y',k);zeros(k,1)]; 6 end 7 yi=Y(1,1); 8 for i=1:n-1 9 z=t; 10 for k=1:i-1 11 z=z*(t-k); 12 end 13 yi=yi+Y(1,i+1)*z/prod([1:i]); 14 end
運行實例：
等距節點的牛頓向前插值代碼：
1 function yi=newtonint2(x,y,xi) 2 n=length(x);h=x(n)-x(n-1);t=(x(n)-xi)/h; 3 n=length(y);Y=zeros(n);Y(:,1)=y'; 4 for k=1:n-1 5 Y(:,k+1)=[zeros(k,1);diff(y',k)]; 6 end 7 h=x(n)-x(n-1);t=(x(n)-xi)/h;yi=Y(n,1); 8 for i=1:n-1 9 z=t; 10 for k=1:i-1 11 z=z*(t-k); 12 end 13 yi=yi+Y(n,i+1)*(-1)^i*z/prod([1:i]); 14 end
運行示例：
3.使用4次牛頓插值多項式插值，並作圖：
解：由4次牛頓插值多項式，
求上述多項式的系數：（修改newtonint.m代碼，得到差商表），代碼如下：
1 function B=newtonint4(x,y) 2 m=length(x);n=length(y); 3 if m~=n 4 error('向量x與y的長度必須一致'); 5 end 6 A=zeros(n); 7 A(:,1)=y; 8 for j=2:n%j為列標 9 for i=1:(n-j+1) %i為行標 10 A(i,j)=(A(i+1,j-1)-A(i,j-1))/(x(i+j-1)-x(i));%計算差商表 11 end 12 end 13 B=A; 14 end
代入數據得到差商表：
| 0.98 |
-0.3 |
-0.625 |
-0.2083 |
-0.5208 |
| 0.92 |
-0.55 |
-0.75 |
-0.625 |
0 |
| 0.81 |
-0.85 |
-1.125 |
0 |
0 |
| 0.64 |
-1.3 |
0 |
0 |
0 |
| 0.38 |
0 |
0 |
0 |
0 |
已知，第一行的便是插值多項式的系數，代入插值多項式：
並作出圖像：
1 x0=[0.2 0.4 0.6 0.8 1.0]; 2 y0=[0.98 0.92 0.81 0.64 0.38]; 3 plot(x0,y0,'b-o') 4 hold on 5 k=0:1:10; 6 x=0.2+0.08*k; 7 for i=1:1:11 8 y(i)=0.98-0.3*(x(i)-0.2)-0.625*(x(i)-0.2)*(x(i)-0.4)-0.2083333*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.520833333*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8); 9 end 10 plot(x,y,'r-o'); 11 legend('原圖像','4次插值圖像');
小結：
在編寫牛頓插值的代碼時，我遇到了超出元組索引的問題。我在MATLAB的提示下（它的提示是英語），如圖：
這個f是使用迭代來求差商的，但是出現了問題。我根據它的提示創建了一個全零數組用於存儲運算得到的差商，在某種程度上解決了這個問題。
在解決第3題時，我特意編寫了一個算差商的程序和一個4次牛頓插值多項式代入數據畫圖的程序。差商的程序是修改第2題的牛頓插值程序得到的，這在一定程度上說明，一個程序的功能是可以分開的同時也可以寫在一起的。但在寫4次多項式代入畫圖的程序時，並沒有參考的我，只能回看書本關於4次牛頓插值的知識，我得到了這個牛頓插值多項式的公式，並發現它的關鍵就是每一項的系數，而那些系數就是算得的差商，所以，很快，我就寫出了4次多項式代入畫圖的程序。很開心的是，算得的多項式擬合得很好。