實驗目的
(1)熟悉拉格朗日插值法、分段線性插值、三次樣條插值等多項式的應用,了解各方法的使用范圍。
(2)熟練掌握三次樣條插值不同條件下的使用和多項式插值的高次插值的使用。
實驗要求
實驗步驟要有模型建立,模型求解、結果分析。
實驗要求
(1)在區間[-1.,1]上分別取n=10,n=20用兩組等距節點對龍格函數f(x)=1/(1+25x2)做拉格朗日插值及三次樣條插值(第一邊界條件,端點的一階導數值),對每個n值,分別畫出插值函數及f(x)的圖像。
(2)已知數據
x |
0 |
1 |
4 |
9 |
16 |
25 |
36 |
49 |
64 |
y |
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
可以得到平方根函數的近似,在區間[0,64]上作圖。
(1)用這9個點作8次拉格朗日多項式插值L8(x).
(2)用三次樣條插值(第一邊界條件)程序求S(x)。S'(x0)=f0',S'(xn)=fn',從得到結果看在【0,64】上,哪個插值更精確;在區間【0,1】上,n=8兩種插值哪個更精確?
實驗步驟
(1)解:本報告編寫拉格朗日插值法代碼得到拉格朗日每個n值的情形,使用MATLAB自帶的interp1()函數實現三樣條插值。拉格朗日代碼如下,

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
n=10時代碼如下,(n=20只需改第2行的值)

1 %n=10 2 n=10; 3 x=linspace(-1,1,n); 4 y=1./(1+25*x.^2); 5 x1=linspace(-1,1,100); 6 y1=1./(1+25*x1.^2); 7 %調用插值函數 8 p=Lagrange(x,y,x1);%拉格朗日 9 p1=interp1(x,y,x1,'spline');%三樣條插值 10 figure,plot(x1,y1,'m','LineWidth',2); 11 hold on 12 plot(x1,p,'g','LineWidth',2); 13 hold on 14 plot(x1,p1,'b','LineWidth',2); 15 legend('龍格函數','拉格朗日插值函數','三樣條插值函數'); 16 %grid on; 17 title('n=10的插值函數及原函數圖像'); 18 xlabel('x軸'); 19 ylabel('y軸'); 20
運行程序,結果如下圖
當n=10時,可見三樣條插值較拉格朗日插值效果更好。
當n=20時,運行程序,(為了區分圖像,本報告對龍格函數線條做了修改)
由上圖可見,對龍格函數插值使用三樣條插值較拉格朗日插值效果更好。
(2)解:本報告編寫相應的拉格朗日代碼和三樣條插值代碼,並畫出[0,64]和[0,1]的這三者的函數圖像及誤差圖。拉格朗日代碼

1 function [yy,p]=lagrange1(x1,y1,xx) 2 %本程序為Lagrange1插值,其中x1,y1 3 %為插值節點和節點上的函數值,輸出為插值點xx的函數值, 4 %xx可以是向量。 5 n=length(x1); 6 m=length(y1); 7 if m~=n 8 error('輸入有誤!!'); 9 end 10 if nargin<3 11 xx=linspace(min(x1),max(x1),n*10); 12 syms x 13 for i=1:n 14 t=x1; 15 t(i)=[]; 16 L(i)=prod((x-t)./(x1(i)-t));% L向量用來存放插值基函數 17 end 18 u=sum(L.*y1); 19 p=simplify(u); % p是簡化后的Lagrange插值函數(字符串) 20 yy=subs(p,x,xx); %p是以x為自變量的函數,並求xx處的函數值 21 end
求三樣條插值S(x)代碼

1 %三樣條擬合函數 2 %輸入:x1,y1插值點 3 %輸出:方程與函數值 4 function [yy,s3]=sanci(x1,y1) 5 syms x; 6 n=length(x1); 7 m=length(y1); 8 if n~=m 9 error('輸入參數有誤'); 10 end 11 p=polyfit(x1,y1,3); 12 s3=p(1)+p(2).*x+p(3).*x.^2+p(4).*x.^3; 13 x2=linspace(min(x1),max(x1),x1(n)); 14 yy=polyval(p,x2); 15 %plot(x2,yy); 16 end
解題1及題2的部分,解題1部分圖
L8(x)= -(x*(143*x^7 - 29260*x^6 + 2366546*x^5 - 97191380*x^4 + 2171047879*x^3 -26340674360*x^2 + 166253376432*x - 577880352000))/435891456000。
解題2部分截圖
S(x)= (3500*x^3)/6927 + (2762705724454173*x^2)/9007199254740992 - (3397610023959411*x)/576460752303423488 + 6795786867330573/147573952589676412928
作[0,64]的函數圖像,代碼

1 %9個點 2 x=[0 1 4 9 16 25 36 49 64]; 3 y=[0 1 2 3 4 5 6 7 8]; 4 x1=linspace(0,64,64); 5 %調用插值函數 6 [p,S]=lagrange1(x,y,x1);%拉格朗日 7 [p1,s3]=sanci(x,y);%三樣條插值 8 figure,plot(x1,sqrt(x1),'rs','LineWidth',2); 9 hold on 10 plot(x1,p,'g','LineWidth',2); 11 hold on 12 plot(x1,p1,'b','LineWidth',2); 13 %設置圖例,命名並設置放在左上角位置 14 legend('平方根函數','拉格朗日插值函數','三樣條插值函數','Location','northwest'); 15 grid on; 16 title('[0,64]的n=9的插值函數及原函數圖像'); 17 xlabel('x軸'); 18 ylabel('y軸'); 19
運行示例:
在[0,64]的插值中,顯然三樣條插值更好。那么在[0,1],n=8條件下的插值呢,對上述程序稍作修改,運行結果如下,
顯然,在[0,1],n=8的條件插值中,拉格朗日的效果更好,代碼如下

1 %8個點 2 x=linspace(0,1,8); 3 y=sqrt(x); 4 x1=linspace(0,1,20); 5 %調用插值函數 6 [p,S]=lagrange1(x,y,x1);%拉格朗日 7 [p1,s3]=sanci(x,y);%三樣條插值 8 figure,plot(x1,sqrt(x1),'rs','LineWidth',2); 9 hold on 10 plot(x1,p,'g','LineWidth',2); 11 hold on 12 plot(x1,p1,'bo','LineWidth',2); 13 %設置圖例,命名並設置放在右下角位置 14 legend('平方根函數','拉格朗日插值函數','三樣條插值函數','Location','southeast'); 15 grid on; 16 17 title('[0,1]的n=8的插值函數及原函數圖像'); 18 xlabel('x軸'); 19 ylabel('y軸');
小結
在解第2題時,發現似乎在小區間的插值,拉格朗日的插值效果要比三次樣條的插值效果更好,到大區間時,三次樣條的插值效果要比拉格朗日的插值效果更好。