数值分析作业姓名王建忠学号132080202006学院能源与动力工程专业机械电子工程2013年12月16日grange插值多项式程序function f=nalagr(x,y,xx)%x为节点向量%y为节点函数值%xx是插值点syms sif(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等!');return;endf=0.0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(s-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(s-x(j))/(x(i)-x(j));%计算拉格朗日基函数end;f=f+l;%计算拉格朗日插值函数simplify(f);if(i==n)if(nargin==3)f=subs(f,'s');%计算插值点的函数值elsef=collect(f);%将插值多项式展开f=vpa(f,6);%将插值多项式的系数化成6位精度的小数endendend>>x=[-2,-1,0,1];%已知节点向量y=[3,1,1,6];%节点函数值向量f=nalagr(x,y)f=0.5*s^3+ 2.5*s^2+ 2.0*s+ 1.0>>f=nalagr(x,y,0)f=1>>2.牛顿插值多项式程序function[p2,z]=newTon(x,y,t)%输入参数中x,y为元素个数相等的向量，t为待估计的点，可以为数字或向量。

%输出参数中p2为所求得的牛顿插值多项式，z为利用多项式所得的t的函数值。

n=length(x);chaS(1)=y(1);for i=2:nx1=x;y1=y;x1(i+1:n)=[];y1(i+1:n)=[];n1=length(x1);s1=0;for j=1:n1t1=1;for k=1:n1if k==jcontinue;elset1=t1*(x1(j)-x1(k));endends1=s1+y1(j)/t1;endchaS(i)=s1;endb(1,:)=[zeros(1,n-1)chaS(1)];cl=cell(1,n-1);for i=2:nu1=1;for j=1:i-1u1=conv(u1,[1-x(j)]);cl{i-1}=u1;endcl{i-1}=chaS(i)*cl{i-1};b(i,:)=[zeros(1,n-i),cl{i-1}];endp2=b(1,:);for j=2:np2=p2+b(j,:);endif length(t)==1rm=0;for i=1:nrm=rm+p2(i)*t^(n-i);endz=rm;elsek1=length(t);rm=zeros(1,k1);for j=1:k1for i=1:nrm(j)=rm(j)+p2(i)*t(j)^(n-i);endz=rm;endend>>x=[-1,1,2,5];y=[-7,7,-4,35];t=1;[u,v]=newTon(x,y,t)u=2-10510%多项式系数v=7%插值点函数值>>3.牛顿迭代法求多项式根%主程序：function[k,x,wuca,yx]=newtondiedai(x0,e)k=1;yx1=fun(x0);yx2=fun1(x0);x1=x0-yx1/yx2;while abs(x1-x0)>ex0=x1;yx1=fun(x0);yx2=fun1(x0);k=k+1;x1=x1-yx1/yx2;endk;x=x1;wuca=abs(x1-x0)/2;yx=fun(x);end%分程序1：function y1=fun(x)y1=x-cos(x);%待求解多项式end%分程序2：function y2=fun1(x)y2=1+sin(x);%函数fun(x)的导数end%分程序2：function y2=fun1(x)y2=1+sin(x);%函数fun(x)的导数end>>[k,x,wuca,yx]=newton(1,10^-5)k=4x=0.7391wuca=3.4480e-10yx=04.三次样条插值函数Matlab程序function[]=spline3(X,Y,dY,x0,m)N=size(X,2);s0=dY(1);sN=dY(2);interval=0.025;disp('x0为插值点')x0h=zeros(1,N-1);for i=1:N-1h(1,i)=X(i+1)-X(i);endd(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);d(N,1)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);for i=2:N-1d(i,1)=6*((Y(1,i+1)-Y(1,i)))/(h(1,i)-(Y(1,i)-Y(1,i-1)))/h(1,i-1) /(h(1,i)+h(1,i-1));endmu=zeros(1,N-1);md=zeros(1,N-1);md(1,N-1)=1;mu(1,1)=1;for i=1:N-2u=h(1,i+1)/(h(1,i)+h(1,i+1));mu(1,i+1)=u;md(1,i)=1-u;endp(1,1)=2;q(1,1)=mu(1,1)/2;for i=2:N-1p(1,i)=2-md(1,i-1)*q(1,i-1);q(1,i)=mu(1,i)/p(1,i);endp(1,N)=2-md(1,N-1)*q(1,N-1);y=zeros(1,N);y(1,1)=d(1)/2;for i=2:N y(1,i)=(d(i)-md(1,i-1)*y(1,i-1))/p(1,i);endx=zeros(1,N);x(1,N)=y(1,N);for i=N-1:-1:1x(1,i)=y(1,i)-q(1,i)*x(1,i+1);endfprintf('M为三对角方程的解\n');M=xfprintf('\n');syms t;digits(m);for i=1:N-1pp(i)=M(i)*(X(i+1)-t)^3/(6*h(i))+M(i+1)*(t-X(i))^3/(6*h(i))+(Y(i)-M(i) *h(i)^2/6)*(X(i+1)-t)/h(i)+(Y(i+1)-M(i+1)*h(i)^2/6)*(t-X(i))/h(i);pp(i)=simplify(pp(i));coeff=sym2poly(pp(i));if length(coeff)~=4tt=coeff(1:3);coeff(1:4)=0;coeff(2:4)=tt;endif x0>X(i)&x0<X(i+1)L=i;y0=coeff(1)*x0^3+coeff(2)*x0^2+coeff(3)*x0+coeff(4);endval=X(i):interval:X(i+1);for k=1:length(val)fval(k)=coeff(1)*val(k)^3+coeff(2)*val(k)^2+coeff(3)*val(k)+coeff(4); endif mod(i,2)==1plot(val,fval,'r+')else plot(val,fval,'b.')endhold onclear val fvalans=sym(coeff,'d');ans=poly2sym(ans,'t');fprintf('在区间[%f,%f]内\n',X(i),X(i+1));fprintf('三次样条函数S(%d)=',i);pretty(ans);endfprintf('x0所在区间为[%f,%f]\n',X(L),X(L+1));fprintf('函数在插值点x0=%f的值为\n',x0);y0>>X=[12346];Y=[00.6931471.098611.38629 1.79176];dY=[11/6];x0=5;m=5;spline3(X,Y,dY,x0,m)x0为插值点x0=5M为三对角方程的解M=-2.1818 2.52260.01970.3019-0.2050在区间[1.000000,2.000000]内三次样条函数S(1)=320.78407t- 3.4431t+ 5.5341t- 2.875在区间[2.000000,3.000000]内三次样条函数S(2)=32-0.41715t+ 3.7642t-10.49t+9.9528在区间[3.000000,4.000000]内三次样条函数S(3)=320.047026t-0.41338t+ 1.4414t-0.77484在区间[4.000000,6.000000]内三次样条函数S(4)=32-0.042241t+0.65781t- 3.1651t+ 6.2252 x0所在区间为[4.000000,6.000000]函数在插值点x0=5.000000的值为y0=1.5648。