MATLAB拉格郎日插值法与牛顿插值法构造插值多项式

作业报告使用MATLAB编程2.1Newton插值：clearclcx=[0.2 0.4 0.6 0.8 1];n=length(x);y=zeros(n);y(:,1)=[0.98 0.92 0.81 0.64 0.38]; a=x(1);b=x(end);for j=1:n-1for k=j:n-1temp=y(k+1,j)-y(k,j);y(k+1,j+1)=temp/(x(k+1)-x(k+1-j)); endc(j)=y(j,j);c(j+1)=y(j+1,j+1);endp=c(1);q=1;syms Xfor i=2:nq=q*(X-x(i-1));p=p+c(i)*q;endp=simple(p)for i=1:301t(i)=a+(b-a)/300*(i-1);Nn(i)=eval(subs(p,'X','t(i)')); endplot(t,Nn,'b')hold ongrid onlegend('Newton')title('Newton')xlabel('x')ylabel('f(x)')图像：三次样条插值：M文件spline33.mfunction [yy b c d]=spline3(x,y,xx,flag,vl,vr)if length(x)~=length(y)error('输入数据应成对');endn=length(x);a=zeros(n-1,1);b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,1);fori=1:n-1a(i)=y(i);dx(i)=x(i+1)-x(i);dy(i)=y(i+1)-y(i);endfori=2:n-1A(i,i-1)=dx(i-1);A(i,i)=2*(dx(i-1)+dx(i));A(i,i+1)=dx(i);B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1)); endif flag==0A(1,1)=1;A(n,n)=1;end% --------------------------------- %% 端点一阶导数条件%if flag==1A(1,1)=2*dx(1);A(1,2)=dx(1);A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1);B(1,1)=3*(dy(1)/dx(1)-vl);B(n,1)=3*(vr-dy(n-1)/dx(n-1));end% --------------- %% 端点二阶导数条件%if flag==2A(1,1)=2;A(n,n)=2;B(1,1)=vl;B(n,1)=vr;end% --------------- %c=A\B;fori=1:n-1d(i)=(c(i+1)-c(i))/(3*dx(i));b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3; end[mmnn]=size(xx);yy=zeros(mm,nn);fori=1:mm*nnfor ii=1:n-1if xx(i)>=x(ii)&xx(i)<x(ii+1)j=ii;break;elseif xx(i)==x(n)j=n-1;endendyy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3; endend主程序：clear;clc;format short g;x1=[0.2 0.4 0.6 0.8 1.0]';y1=[0.98 0.92 0.81 0.64 0.38]';u1=0;un=0;xx1=[x1(1):0.001:x1(end)]';[yy1 b1 c1 d1]=spline33(x1,y1,xx1,0,u1,un);fprintf('\t\tb1\t\tc1\t\td1\n');disp([b1 c1(1:end-1,1) d1]);plot(x1,y1,'bo',xx1,yy1,'--k');legend('spline')title('spline')grid on;2.3（1）8次拉格朗日插值M文件language1.mfunction f=language(x,y,x0) syms t l;if (length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等！') return;endh=sym(0);fori=1:nl=sym(y(i));for j=1:i-1l=l*(t-x(j))/(x(i)-x(j)); end;for j=i+1:nl=l*(t-x(j))/(x(i)-x(j)); end;h=h+l;endsimplify(h);if(nargin==3)f=subs(h,t,x0);elsef=collect(h);f=vpa(f,6);end主程序：x=[0 1 4 9 16 25 36 49 64];y=[0 1 2 3 4 5 6 7 8];f=language1(x,y)plot(x,y)legend('language')title('language')grid on;（2）三次样条差值：M文件spline33.mfunction [yy b c d]=spline3(x,y,xx,flag,vl,vr) if length(x)~=length(y)error('输入数据应成对');endn=length(x);a=zeros(n-1,1);b=a;d=a;dx=a;dy=a;A=zeros(n);B=zeros(n,1);fori=1:n-1a(i)=y(i);dx(i)=x(i+1)-x(i);dy(i)=y(i+1)-y(i);endfori=2:n-1A(i,i-1)=dx(i-1);A(i,i)=2*(dx(i-1)+dx(i));A(i,i+1)=dx(i);B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));endif flag==0A(1,1)=1;A(n,n)=1;end% --------------------------------- %% 端点一阶导数条件%if flag==1A(1,1)=2*dx(1);A(1,2)=dx(1);A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1);B(1,1)=3*(dy(1)/dx(1)-vl);B(n,1)=3*(vr-dy(n-1)/dx(n-1));end% --------------- %% 端点二阶导数条件%if flag==2A(1,1)=2;A(n,n)=2;B(1,1)=vl;B(n,1)=vr;end% --------------- %c=A\B;fori=1:n-1d(i)=(c(i+1)-c(i))/(3*dx(i));b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;end[mmnn]=size(xx);yy=zeros(mm,nn);fori=1:mm*nnfor ii=1:n-1if xx(i)>=x(ii)&xx(i)<x(ii+1)j=ii;break;elseif xx(i)==x(n)j=n-1;endendyy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3;endend主程序：clear;clc;format short g;x1=[0 1 4 9 16 25 36 49 64 ]';y1=[0 1 2 3 4 5 6 7 8]';u1=0;un=0;xx1=[x1(1):0.001:x1(end)]';[yy1 b1 c1 d1]=spline33(x1,y1,xx1,1,u1,un); fprintf('\t\tb1\t\tc1\t\td1\n');disp([b1 c1(1:end-1,1) d1]);plot(x1,y1,'bo',xx1,yy1,'--k');legend('spline')title('spline')grid on;4-1(1)M文件function.m:Function y=fun(x);y=sqrt(x)*log(x);主程序Clear；clc；h=0.001；n=1/h；t=0；s1=0；s2=0；fori=1:n-1t=t+function(i*h);endT=h/2*(0+2*t+f(1));T=vpa(T,7)For i=0:n-1s1=s1+function(h/2+i*h);endfori=1:n-1s2=s2+function(i*h);endS=h/6*(0+4*s1+2*s2+f(1));S=vpa(S,7) 综上，得解。

Matlab插值法实验⽬的：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);3if m~=n4 error('向量x与y的长度必须⼀致');5 end6 s=0;7for k=1:n8 t=ones(1,p);9for j=1:n10if j~=k11 t=t.*(xi-x(j))./(x(k)-x(j));12 end13 end14 s=s+t.*y(k);15 end16 yi=s;17 endLagrange 运⾏⽰例： ⽜顿插值法代码：1 function yi=newtonint(x,y,xi)2 m=length(x);n=length(y);3if m~=n4 error('向量x与y的长度必须⼀致');5 end6 A=zeros(n);7 A(:,1)=y;8for j=2:n%j为列标9for 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 end12 end13 %根据差商表,求对应的⽜顿插值多项式在x=xi处的值yi14 N(1)=A(1,1);15for j=2:n16 T=1;17for i=1:j-118 T=T*(xi-x(i));19 end20 N(j)=A(1,j)*T;21 end22 yi=sum(N); %将x=xi带⼊⽜顿插值多项式，得到的yi的值23 %A 输出差商表24 endnewtonint 运⾏实例： 等距节点的⽜顿向后插值代码：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';4for k=1:n-15 Y(:,k+1)=[diff(y',k);zeros(k,1)];6 end7 yi=Y(1,1);8for i=1:n-19 z=t;10for k=1:i-111 z=z*(t-k);12 end13 yi=yi+Y(1,i+1)*z/prod([1:i]);14 endnewtonint1 运⾏实例： 等距节点的⽜顿向前插值代码：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';4for k=1:n-15 Y(:,k+1)=[zeros(k,1);diff(y',k)];6 end7 h=x(n)-x(n-1);t=(x(n)-xi)/h;yi=Y(n,1);8for i=1:n-19 z=t;10for k=1:i-111 z=z*(t-k);12 end13 yi=yi+Y(n,i+1)*(-1)^i*z/prod([1:i]);14 endnewtonint2 运⾏⽰例： 3.使⽤4次⽜顿插值多项式插值，并作图： 解：由4次⽜顿插值多项式， 求上述多项式的系数：（修改newtonint.m代码，得到差商表），代码如下：1 function B=newtonint4(x,y)2 m=length(x);n=length(y);3if m~=n4 error('向量x与y的长度必须⼀致');5 end6 A=zeros(n);7 A(:,1)=y;8for j=2:n%j为列标9for 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 end12 end13 B=A;14 endnewtonint4 代⼊数据得到差商表：0.98-0.3-0.625-0.2083-0.52080.92-0.55-0.75-0.62500.81-0.85-1.125000.64-1.30000.380000 已知，第⼀⾏的便是插值多项式的系数，代⼊插值多项式： 并作出图像：1 x0=[0.20.40.60.81.0];2 y0=[0.980.920.810.640.38];3 plot(x0,y0,'b-o')4 hold on5 k=0:1:10;6 x=0.2+0.08*k;7for i=1:1:118 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 end10 plot(x,y,'r-o');11 legend('原图像','4次插值图像');plot3⼩结： 在编写⽜顿插值的代码时，我遇到了超出元组索引的问题。

拉格朗⽇插值和⽜顿插值matlab1. 已知函数在下列各点的值为0.20.40.60.8 1.00.980.920.810.640.38⽤插值法对数据进⾏拟合，要求给出Lagrange插值多项式和Newton插值多项式的表达式，并计算插值多项式在点的值。

程序：x=[0.2 0.4 0.6 0.8 1.0];y=[0.98 0.92 0.81 0.64 0.38];x0=[0.2 0.28 0.44 0.76 1 1.08];[f,f0]=Lagrange(x,y,x0)function [f,f0] = Lagrange(x,y,x0)%求已知数据点的Lagrange插值多项式f，并计算插值多项式f在数据点x0的函数值f0syms t;n = length(x);f = 0.0;for i = 1:nl = y(i);for j = 1:i-1l = l*(t-x(j))/(x(i)-x(j));end;for j = i+1:nl = l*(t-x(j))/(x(i)-x(j));end;f = f + l;simplify(f);if(i==n)f0 = subs(f,'t',x0);f = collect(f);f = vpa(f,6);endend结果：>> Untitled3f =- 0.520833*t^4 + 0.833333*t^3 - 1.10417*t^2 + 0.191667*t + 0.98f0 =[ 49/50, 60137/62500, 56377/62500, 42497/62500, 19/50, 15017/62500]⽜顿：%y为对应x的值，A为差商表，C为多项式系数，L为多项式%X为给定节点，Y为节点值，x为待求节点function[y,A,C,L] = newton(X,Y,x,M)n = length(X);m = length(x);for t = 1 : mz = x(t);A = zeros(n,n);A(:,1) = Y';s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;for j = 2 : nfor i = j : nA(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));endq1 = abs(q1*(z-X(j-1)));c1 = c1 * j;endC = A(n, n); q1 = abs(q1*(z-X(n)));for k = (n-1):-1:1C = conv(C, poly(X(k)));d = length(C);C(d) = C(d) + A(k,k);endy(t) = polyval(C,z);endL = poly2sym(C);x=[0.2 0.4 0.6 0.8 1.0];y=[0.98 0.92 0.81 0.64 0.38];x0=[0.2 0.28 0.44 0.76 1 1.08];m=1;[y,A,C,L]=newton(x,y,x0,m)结果：y =0.9800 0.9622 0.9020 0.6800 0.3800 0.2403A =0.9800 0 0 0 00.9200 -0.3000 0 0 00.8100 -0.5500 -0.6250 0 00.6400 -0.8500 -0.7500 -0.2083 00.3800 -1.3000 -1.1250 -0.6250 -0.5208C =-0.5208 0.8333 -1.1042 0.1917 0.9800L =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/502. 在区间上分别取，⽤两组等距节点对Runge函数作多项式插值（Lagrange插值和Newton插值均可），要求对每个值，分别画出插值多项式和函数的曲线。

大家都知道插值在数学建模中很重要，现在介绍几种常用插值下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite插值和三次样条插值。

1.拉格朗日多项式插值拉格朗日插值就是给定n个数，让你用不超过n-1次的多项式你逼近它，当然这n个点要能满足多项式。

这是一种最基本的思想，计算很简单，先计算n个基函数，基函数可以自己上网搜一下，因为这里打出公式有点麻烦。

然后就是把每个点的y值乘以他的基函数，把这n个式子相加，最后化简就ok了。

下面我把代码写出来，我这些代码全是自己写的，注释比较详细，这里只以lagrange为例，其余都放在附件里了。

%定义myLagrange函数，参数为向量x，y，由用户调用该函数时输入function L=myLagrange (x,y)%n 插值结点的个数n=length(x);%L myLagrange函数计算的多项式系数行列式L=zeros(1,n);%%使用双重for循环，第一个for循环是fori=1:n%aa=1;%ww=1;%for循环for j=1:n%如果i不等于jif j~=i%累加法计算aa=a*(x(i)-x(j));%用向量乘法函数conv计算ww=conv(w,[1,-x(j)]);%if语句结束符end%第二个for循环结束符end%递归法计算L，其中y(i)/a*w表示第i个元素L=y(i)/a*w+L;%第一个for结束符end没错，就这么几句代码，所以很简单的。

2.牛顿插值牛顿插值其实是为了解决拉格朗日插值不能增加新的点来说的。

拉格朗日插值只能接受给定的那么多点，了然后插值。

如果你想再加一个点，它会重新开始计算，这个很费时间和内存。

因此牛顿插值就诞生了。

了解牛顿插值前要学习下差商和差分两个简单的概念。

Newton 插值的优点是：每增加一个节点，插值多项式只增加一项，即因而便于递推运算。

而且Newton 插值的计算量小于Lagrange 插值。