搜索
数值分析计算实习题第二章2-1程序：clear;clc;x1=[0.20.40.60.81.0];y1=[0.980.920.810.640.38];n=length(y1);c=y1(:);forj=2:n%求差商fori=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsymsxdfd;df(1)=1;d(1)=y1(1);fori=2:n%求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5)%P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1,'variational');%调用三次样条函数q=pp.coefs;disp('三次样条函数');fori=1:4S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];S=vpa(collect(S),5)endx2=0.2:0.08:1.08;dot=[121112];figureezplot(P4,[0.2,1.08]);holdony2=fnval(pp,x2);x=x2(dot);y3=eval(P4);y4=fnval(pp,x2(dot));plot(x2,y2,'r',x2(dot),y3,'b*',x2(dot),y4,'co');title('4次牛顿插值及三次样条');结果如下：4次牛顿插值多项式P4=-0.52083*x^4+0.83333*x^3-1.1042*x^2+0.19167*x+0.98三次样条函数x∈[0.2,0.4]时，S=-1.3393*x^3+0.80357*x^2-0.40714*x+1.04x∈[0.4,0.6]时，S=0.44643*x^3-1.3393*x^2+0.45*x+0.92571x∈[0.6,0.8]时，S=-1.6964*x^3+2.5179*x^2-1.8643*x+1.3886x∈[0.8,1.0]时，S=2.5893*x^3-7.7679*x^2+6.3643*x-0.80571输出图如下2-3（1）程序：clear;clc;x1=[01491625364964];y1=[012345678];%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;fori=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);fori=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:));%wn+1/(x-xk)end三次样条插值曲线4次牛顿插值曲线Dw=zeros(1,n);fori=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));%系数endp=Dw*q;symsxL8;fori=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('8次拉格朗日插值');L8=vpa(collect((sum(L8))),5)xi=0:64;yi=polyval(p,xi);figureplot(xi,yi,x1,y1,'r*');holdontitle('8次拉格朗日插值');结果如下：8次拉格朗日插值L8=-3.2806e-10*x^8+6.7127e-8*x^7-5.4292e-6*x^6+0.00022297*x^5-0.0049807*x^4+0.060429*x^3-0.38141*x^2+1.3257*x输出图如下：第五章4-1（3）程序：clc;clear;y=@(x)sqrt(x).*log(x);a=0;b=1;tol=1e-4;p=quad(y,a,b,tol);fprintf('采用自适应辛普森积分结果为：%d\n',p);结果如下：采用自适应辛普森积分结果为：-4.439756e-01第九章9-1(a)程序：clc;clear;a=1;b=2;%定义域h=0.05;%步长n=(b-a)/h;y0=1;%初值f=@(x,y)1/x^2-y/x;%微分函数Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解fori=1:nxn=Xn(i);xnn=Xn(i+1);yn=Yn(i);yp=yn+h*f(xn,yn);yc=yn+h*f(xnn,yp);yn=(yp+yc)/2;Yn(i+1)=yn;endXn=Yn;%以下根据经典四阶R-K法公式求解fori=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp('改进欧拉法四阶经典R-K法');disp([Xn'Yn'])结果如下：改进欧拉法四阶经典R-K法110.998870.998850.995770.99780.991140.996940.985320.996340.978570.996030.971110.996060.963110.996450.95470.997230.945980.998410.9370510.927981.0020.918831.00440.909641.00730.900451.01060.891291.01430.882181.01840.873151.02290.864211.02780.855381.03310.846651.0388（b）程序：clc;clear;a=0;b=1;%定义域H=[0.10.0250.01];%步长y0=1/3;%初值f=@(x,y)-50*y+50*x^2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.^2;%准确解Ym=zeros(1,11);forj=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值fori=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;endfork=1:11m=0.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为：%d\n',h);disp('四阶经典R-K法准确解误差');disp([Ym'Y'delta'])end结果如下：步长为：1.000000e-01四阶经典R-K法准确解误差0.333330.3333304.60550.0122464.593263.0620.04001563.022864.050.09863.96118440.16118431.6235e+050.251.6235e+052.2256e+060.362.2256e+063.0509e+070.493.0509e+074.1823e+080.644.1823e+085.7333e+090.815.7333e+097.8594e+1017.8594e+10步长为：2.500000e-02四阶经典R-K法准确解误差0.333330.3333300.0130150.0122460.000768940.0400630.0400154.82e-050.0900370.093.6857e-050.160040.163.6723e-050.250040.253.6722e-050.360040.363.6722e-050.490040.493.6722e-050.640040.643.6722e-050.810040.813.6722e-05113.6722e-05步长为：1.000000e-02四阶经典R-K法准确解误差0.333330.3333300.0122560.0122469.5673e-060.0400160.0400157.8252e-070.0900010.096.6347e-070.160.166.6226e-070.250.256.6225e-070.360.366.6225e-070.490.496.6225e-070.640.646.6225e-070.810.816.6225e-07116.6225e-07由结果可知，步长越小，结果越精确。9-3主程序：clc;clear;options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-41e-4]);[X,Y]=ode45(@rigid,[00.004],[110],options);subplot(1,2,1);plot(X,Y(:,1),'-',X,Y(:,2),'-.',X,Y(:,3),'.');title('四阶R-K方法');legend('y1','y2','y3');[X,YY]=ode23t(@rigid,[00.004],[110],options);subplot(1,2,2);plot(X,YY(:,1),'-',X,YY(:,2),'-.',X,YY(:,3),'.');title('梯形法');legend('y1','y2','y3');rigid函数functiondy=rigid(~,y)dy=zeros(3,1);%acolumnvectordy(1)=-0.013*y(1)-1000*y(1)*y(2);dy(2)=-2500*y(2)*y(3);dy(3)=-0.013*y(1)-1000*y(1)*y(2)-2500*y(2)*y(3);结果如图：