2-MATLAB在数值分析中的应用

-1-第3章MATLAB在数值分析中的应用-2-t=dat(:,1);CA1=dat(:,2);CA2=dat(:,3);plot(t,CA1,'o',t,CA2,'^')axis([0t(end)+0.5min(min(dat(:,2:3)))-0.1max(max(dat(:,2:3)))+0.05])xlabel('Time(min)')ylabel('C_A_1,C_A_2(kmol/m^3)')legend('C_A_1','C_A_2')title('Concentrationprofiles')问题的提出01234567890.20.250.30.350.40.450.5-3-在平面上给定n个点(xk,yk),可以唯一确定一个最多n-1次的多项式通过这些点，这个多项式叫插值多项式插值和拟合P(xk)=yk,k=1,2,…,n1()njkkjkkjxxPxyxxLagrange插值多项式-4-插值多项式例子x0=0:3;y0=[-5–6–116];disp([x0;y0]))16()6()2)(1()1()2()3)(1()6()2()3)(2()5()6()3)(2)(1()(xxxxxxxxxxxxxP52)(3xxxPx0=[0123]y0=[-5-6-116]1()njkkjkkjxxPxyxx-5-Polyinterp(Lagrange插值形式)1()njkkjkkjxxPxyxxfunctiony=polyinterp(x0,y0,x)n=length(x0);m=length(x);z=x(i);s=0.0;p=1.0;ifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endforj=1:nends=s+p*y0(k);fork=1:nendy(i)=s;fori=1:mend-6-12121()...nnnnPxcxcxcxc1211111122222212111nnnnnnnnnnnxxxcycyxxxcyxxxLagrange插值多项式基函数最后一列全为1，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积，这种矩阵称为范得蒙(Vandermonde)矩阵。在MATLAB中，可用函数vander(V)生成以向量V为基础向量生成。x0=0:3;y0=[-5–6–116];A=vander(x0);X=A\B;对A*X=BX=A\Bx0=[0123]B=y0'=[-5-6-116]'x=-.25:.01:3.25;y1=polyval(X,x)-7-x0=0:3;y0=[-5-6-116];x=-.25:.001:3.25;%--------------------------------ticA=vander(x0);X=A\y0';y1=polyval(X,x);t1=toc%--------------------------------ticy2=polyinterp(x0,y0,x);t2=toc%--------------------------------plot(x0,y0,'o',x,y1,'-',x,y2,'b')legend('data','vandmethod','polyinterp')-8--0.500.511.522.533.5-10-50510152025datavandmethodpolyinterp-9-x00=data(:,1);y00=data(:,2);x0=da00(:,1);y0=da00(:,2);x=0:.1:9;%--------------------------------ticA=vander(x0);X=A\y0;y1=polyval(X,x);t1=toc%--------------------------------ticy2=polyinterp(x0,y0,x);t2=toc%--------------------------------plot(x00,y00,'*',x0,y0,'o',x,y1,'-',x,y2,'b')legend('da00','data','vandmethod','polyinterp')-10-Polyinterp(symbol)symx=sym(‘x’)P=polyinterp(x0,y0,symx)pretty(P)P=simplify(P)P=x^3-2*x-5simplify(S)simplifieseachelementofthesymbolicmatrixS8/3*z*(x-1)*(x-2)+1/2*x*(x-1)*(x-3)-3*x*(x-2)*(x-3)-5/6*(-x+1)*(x-2)*(x-3)symx=sym(‘x')createsthesymbolicvariablewithname‘x'pretty(S)printsthesymbolicexpressionSinaformatthatresemblestype-setmathematics.-11-012345676810121416182022Polyinterp(另外的例子)x0=1:6;y0=[161821171512];x=.75:.05:6.25;y=polyinterp(x0,y0,x);plot(x0,y0,’o’,x,y,’-’)-12-x0=[-5:1:5];y0=1./(1+x0.^2);x=[-5:0.1:5];y=polyinterp(x0,y0,x);y1=1./(1+x.^2);plot(x0,y0,'r--')holdonplot(x,y1,'b-')Lagrange插值多项式在全区间内不收敛的情况-5-4-3-2-1012345-0.500.511.52-13-Hermite插值设x1,x2,…,xn对应的函数值y1,y2,…,yn，一阶导数值y1’,y2’,…,yn’Hermite插值公式niiiiiiiyyyaxxhxy1'])2)([()(21)(nijjjijixxxxhnijijiixxa11-14-functiony=hermite(x0,y0,y1,x)n=length(x0);m=length(x);yy=0.0;h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;ifj~=iendforj=1:nendh=1.0;a=0.0;yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));fori=1:nendy(k)=yy;fork=1:mend'1()[()(2)]niiiiiiiyxhxxayyy21)(nijjjijixxxxhnijijiixxa11-16-分段插值一维插值：可以分为最近插值、线性插值、三次样条插值，分段三次Hermite插值。y=interp1(x0,y0,x)y=interp1(x0,y0,x,’method’)method=nearst:最近插值linear:线性插值（默认值）spline:分段三次样条插值pchip:分段三次Hermite插值-17-0123456710121416182022分段线性插值x=1:6;y=[161821171512];11()()kkkkkkyyyxyxxxx-18-k(x=x0(j))=j;forj=1:n-1endn=length(x0);k=ones(size(x));s=x–x0(k);y=y0(k)+s.*d(k);%Evaluateinterpolant%x0小于x的最近位值%Findsubintervalindices,x0(k)=xx0(k+1)d=diff(y0)./diff(x0);functiony=piecelin(x0,y0,x)%Firstdivideddifference11()()kkkkkkyyyxyxxxxx0=[1:1:3];x=[1:0.5:3];-19-x0=[-5:1:5];y0=1./(1+x0.^2);x=[-5:0.1:5];y=interp1(x0,y0,x,’linear’);y1=1./(1+x.^2);y2=piecelin(x0,y0,x)plot(x0,y0,‘^')holdonplot(x,y,’o’,x,y1,'r',x,y2,’-’)-5-4-3-2-101234500.10.20.30.40.50.60.70.80.91datainterp1originfudf-20-0123456710121416182022x0=1:6;y0=[161821171512];x=.75:.05:6.25;y=interp1(x0,y0,x,’pchip’);plot(x0,y0,’o’,x,y,’-’);-21-三次样条三次样条也是分段三次插值函数。可以给出光滑的插值曲线（面）。样条函数必须满足一阶、二阶导数连续。23,0,1,2,3()()kkkkkkkkPxPxssxxsxxsxxxk1,0,11,2,32626kkkkkkkkkkkkkhPxPxsysdPxPxPxssh-22-x=1:6;y=[161821171512];u=.75:.05:6.25;v=interp1(x,y,u,’spline’);plot(x,y,’o’,u,v,’-’);0123456710121416182022-23-插值二维插值：可以分为最近插值、线性插值、三次样条插值，立方插值。z=interp2(x0,y0,z0,x,y)z=interp2(x0,y0,z0,x,y,’method’)method=nearst:最近插值linear:线性插值（默认值）spline:三次样条插值cubic:立方插值-24-x00=data(:,1);y00=data(:,3);x0=da00(:,1);y0=da00(:,3);x=0:.1:9;ticA=vander(x0);X=A\y0;y1=polyval(X,x);t(1)=toc;ticy2=polyinterp(x0,y0,x);t(2)=toc;ticy3=interp1(x0,y0,x,'linear');t(3)=toc;ticy4=interp1(x0,y0,x,'pchip');t(4)=toc;ticy5=interp1(x0,y0,x,'spline');t(5)=toc;%--------------------------------tsubplot(1,2,1);plot(x00,y00,'*',x0,y0,'o',x,y1,'-',x,y2,'b',x,y5,'k-')legend('18da','interda','vandmethod','polyinterp','spline')subplot(1,2,2);plot(x00,y00,'*',x0,y0,'o',x,y3,'-',x,y4,'b',x,y5,'k-')legend('18da','interda','linear','hermite','spline')-25-t=0.00130.00210.00300.00580.0092-26-t=linspace(0,5,100);y=1-cos(3*t).*exp(-t);plot(t,y,'b');grid;holdon；plot(t,0.95*ones(size(t)),'r');holdoffit=min(find(y