数值分析第五版计算实习题

word完美格式数值分析计算实习题第二章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*xP+0.83333*x^3-1.1042*x^2+0.19167*x+0.98三次样条函数xW[0.2,0.4]时，S二-1.3393*xP+0.80357*xV-0.40714*x+1.04xW[0.4,0.6]时，S二0.44643*x^3-1.3393*x^2+0.45*x+0.92571xW[0.6,0.8]时，S二-1.6964*x^3+2.5179*x^2-1.8643*x+1.3886程序：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:nendDw=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*xV-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);和rintf('采用自适应辛普森积分结果为：%d\n',p);结果如下：采用自适应辛普森积分结果为：-4.439756e-01第九章9-1(a)程序：clc;clear;a=1;b=2；%定义域h=0.05；%步长n=(b-a)/h;y0=1;%初值仁@(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);word完美格式范文范例精心整理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;%初值仁@(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:3范文范例精心整理word完美格式h);准确解误差');O.33333O.333334.6O55O.O1224663.O62O.O4OO15864.O5O.O911844O.161.6235e+O5O.252.2256e+O6O.363.O5O9e+O7O.494.1823e+O8O.645.7333e+O9O.817.8594e+1O1四阶经典R-K法准确解h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn二zeros(1,n);%结果矩阵Yn(1)=yO;%赋初值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=O.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为：%d\n',disp('四阶经典R-K法disp([Ym'Y'delta'])end结果如下：步长为：1.OOOOOOe-O1误差O4.593263.O22863.96118431.6235e+O52.2256e+O63.O5O9e+O74.1823e+O85.7333e+O97.8594e+1O步长为：2.5OOOOOe-O2四阶经典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);