实验名称 插值法实验目的(1)学习并熟练掌握MA TLAB 语言的编程;(2)通过课程实习能够应用MATLAB 软件来计算函数的插值,了解函数插值方法。
实验原理牛顿差商形式多项式P(x)=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+…+f[x0,x1,x2…xn](x-x0)…(x-xn-1) 牛顿插值多项式的余项 Rn(x)=f[x0,x1,x2…xn]wn+1(x) 实验题目{1}已知函数在下列各点的值为i x 0.2 0.4 0.6 0.8 1.0 ()i f x 0.98 0.92 0.81 0.64 0.38试用4次牛顿插值多项式()4P x 及三次样条函数()Q x (自然边界条件)对数据进行插 值。
用图给出{(,i i x y ),i x =0.2+0.08i ,i=0,1,11,10},()4P x 及()Q x 。
①实验过程x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:);for j=2:n %求差商 for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end endsyms x df d ;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式 df(i)=df(i-1)*(x-x1(i-1)); d(i)=c(i-1)*df(i); endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数 pp=csape(x1,y1, 'variational');%调用三次样条函数 q=pp.coefs;q1=q(1,:)*[(x-.2)^3;(x-.2)^2;(x-.2);1]; q1=vpa(collect(q1),5)q2=q(1,:)*[(x-.4)^3;(x-.4)^2;(x-.4);1]; q2=vpa(collect(q2),5)q3=q(1,:)*[(x-.6)^3;(x-.6)^2;(x-.6);1];q3=vpa(collect(q3),5)q4=q(1,:)*[(x-.8)^3;(x-.8)^2;(x-.8);1];q4=vpa(collect(q4),5)%求解并化简多项式②运行结果P4 =- 0.20833333333212067373096942901611*x^4 - 0.20833333333212067373096942901611*x^3 + 0.15833*x^2 + 0.96833*x + 0.782q1 =- 1.3392857142898719757795333862305*x^3 + 0.80357*x^2 - 0.40714285714057041332125663757324*x + 1.04q2 =- 1.3392857142898719757795333862305*x^3 + 1.6071*x^2 - 0.8892857142927823588252067565918*x + 1.1643q3 =- 1.3392857142898719757795333862305*x^3 + 2.4107*x^2 - 1.6928571428579743951559066772461*x + 1.4171q4 =- 1.3392857142898719757795333862305*x^3 + 3.2143*x^2 - 2.8178571428288705646991729736328*x + 1.8629③问题结果4次牛顿差值多项式()4P x = 0.98*x - 0.3*(x - 0.2)*(x - 0.4) - 0.625*(x - 0.2)*(x - 0.4)*(x - 0.6) - 0.20833*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.784三次样条差值多项式()Q x ={2}2.在区间[-1,1]上分别取n=10.20用两组等距节点对龙格函数f(x)=22511x+)做多项式及三次样条差值,对每个n 值,分别画出差值函数及f(x)的函数. 龙格图函数做多项式function L=lagrange(a,b,a0) syms x n=length(a) L=0.0 for i=1:n l=b(i); for j=1:i-1l=l.*(x-a(j))/(a(i)-a(j)); endfor j=i+1:nl=l.*(x-a(j))/(a(i)-a(j)); endL=L+l;simplify(L);endL=collect(L)L=vpa(L,6)L=subs(L,'x',a0);endclear allsubplot(1,2,1);a=linspace(-1,1,10);b=1./(1+25.*a.^2);L=Lagrange(a,b)b0=subs(L,'x',a);plot(a,b0,'rs')hold onplot(a,b)title('n=10ʱµÄ²åÖµº¯Êý'); subplot(1,2,2);a=linspace(-1,1,20);b=1./(1+25.*a.^2);L=Lagrange(a,b)b0=subs(L,'x',a);plot(a,b0,'rs');hold on;plot(a,b)title('n=20ʱµÄ²åÖµº¯Êý');图像:三次样条插值函数插值程序:subplot(1,2,1);a=linspace(-1,1,5);b=1./(1+25.*a.^2);a0=linspace(-1,1,10);b0=spline(a,b,a0);plot(a,b,'rs',a0,b0)hold onplot(a,b,'G')title('n=10时的插值函数'); subplot(1,2,2);a=linspace(-1,1,5);b=1./(1+25.*a.^2);a0=linspace(-1,1,20);b0=spline(a,b,a0);plot(a,b,'rs',a0,b0)hold onplot(a,b,'G')title('n=20时的插值函数');图像:y{3}下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似,在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上,哪个插值更精确;在区间[0,1]上,两种插值哪个更精确?①实验过程(1)拉格朗日插值多项式,求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式.输出结果为Ls =-(x*(143*x^7 - 29260*x^6 + 2366546*x^5 - 97191380*x^4 + 2171047879*x^3 - 26340674360*x^2 + 166253376432*x - 577880352000))/435891456000(2)三次样条插值,程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %三条拟和函数S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x)输出结果为:S =(1998674665313721*x^3)/2251799813685248+(4552380473376719*x^2)/180143985 09481984-(4798224608945675*x)/1152921504606846976+4576150135846985/1475739525 89676412928(3)在区间[0,64]上,分别对这两种插值和标准函数作图,plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线,黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合,因此在区间[0,64]上三次样条插值更精确。
在[0,1]区间上由上图看不出差别,不妨代入几组数据进行比较,取x4=[0:0.2:1] x4=[0:0.2:1];sqrt(x4) %准确值subs(Ls,'x',x4) %拉格朗日插值spline(x1,y1,x4) %三次样条插值运行结果为ans =0 0.4472 0.6325 0.7746 0.8944 1.0000ans =0 0.2504 0.4730 0.6706 0.8455 1.0000ans =0 0.2429 0.4630 0.6617 0.8403 1.0000 从这几组数值上可以看出在[0,1]区间上,拉格朗日插值更精确。
数据拟合和最佳平方逼近。