matlab程序与数值实验.docVIP

MATLAB程序与数值实验 一、 插值法(课本第二章) 1、拉格朗日插值方法 function s=lag(x,y,t) %拉格朗日插值方法，可以同时对多点插值，t可以为向量 syms p; n=length(x);%读取x向量的维数 s=0; for(k=1:n) la=y(k); %构造基函数 for(j=1:k-1) la=la*(p-x(j))/(x(k)-x(j)); end; for(j=k+1:n) la=la*(p-x(j))/(x(k)-x(j)); end; s=s+la; simplify(s); end if(nargin==2) s=subs(s,p,x); s=collect(s);%展开多项式 s=vpa(s,6); else m=length(t); for i=1:m temp(i)=subs(s,p,t(i)); end s=temp; End %for example % x=[pi/4 pi/6 pi/3 pi/2]; %y=[cos(pi/4) cos(pi/6) cos(pi/3) cos(pi/2)]; %t=[-40*pi/180 40*pi/180 50*pi/180 70*pi/180 170*pi/180]; %yt=lag(x,y,t) 2、牛顿插值法 function s=niudun(x,y,t) %，可以同时对多点插值，即t可以为向量 syms p; s=y(1); xishu=0; dxs=1; n=length(x);%读取x向量的维数 %构造牛顿插值方法 for (i=1:n-1) for(j=i+1:n) xishu(j)=(y(j)-y(i))/(x(j)-x(i)); end temp1(i)=xishu(i+1); dxs=dxs*(p-x(i)); s=s+temp1(i)*dxs; y=xishu; end simplify(s) if (nargin==2) s=subs(s,p,x); s=collect(s); s=vpa(s,4); else %读取要插值点的向量长度 %可以直接对多点插值机算 m=length(t); for i=1:m temp(i)=subs(s,p,t(i)); end %得到的是系列插值点的插值结果，即得到的是向量，赋值给s s=temp; end %for example %1、已知零阶Bessel函数f(x)在若干点处的函数值为：计算x在1.5处的近似值 %x=[1.0 1.3 1.6 1.9 2.2]; %y=[0.7651977 0.6200860 0.4554022 0.2818186 0.1103623]; %yt=niudun(x,y,1.5) 插值中的Runge现象 syms x f=1/(1+x^2); x=-5:5; y=subs(f,x); chazhi=niudun(x,y); v=[-5,5,-0.5,2]; ezplot(chazhi),axis(v),grid hold on t=-5:0.05:5; yt=subs(f,t); plot(t,yt,’:’) 4、Hermite插值 function f=Hermite(x,y,dy,t) %Hermite插值，x为插值节点，y为插值节点的函数值，dy为插值节点的一阶导数值，t为被插数据，可以为向量 n=length(x); m=length(t); for k=1:m g(k)=0.0; for i=1:n la=1; lp=0.0; for j=1:n if (j~=i) la=la*(t(k)-x(j))/(x(i)-x(j)); lp=lp+1/(x(i)-x(j)); end end temp1=1-2*(t(k)-x(i))*lp; temp2=y(i)*temp1*la^2; temp3=dy(i)*(t(k)-x(i))*la^2; g(k)=g(k)+temp2+temp3; end end f=g; %fo