1、数值分析计算实习题第二章2-1程序：clear;clc;x1=0.2 0.4 0.6 0.8 1.0;y1=0.98 0.92 0.81 0.64 0.38;n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1)/(x1(i)-x1(i-j+1);endendsyms x df d;df(1)=1;d(1)=y1(1);for i=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
2、) %P4即为 4次牛顿插值多项式,并保留小数点后 5位数pp=csape(x1,y1, variational);%调用三次样条函数q=pp.coefs;disp(三次样条函数);for i=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=1 2 11 12;figureezplot(P4,0.2,1.08);hold ony2=fnval(pp,x2);x=x2(dot);y3=eval(P4);y4=fnval(pp,x2(dot);plot(x2,y2,r,x2(
3、dot),y3,b*,x2(dot),y4,co);title(4次牛顿插值及三次样条);结果如下：4次牛顿插值多项式P4 = - 0.52083*x4 + 0.83333*x3 - 1.1042*x2 + 0.19167*x + 0.98三次样条函数x0.2,0.4时， S = - 1.3393*x3 + 0.80357*x2 - 0.40714*x + 1.04x0.4,0.6时，S = 0.44643*x3 - 1.3393*x2 + 0.45*x + 0.92571x0.6,0.8时，S = - 1.6964*x3 + 2.5179*x2 - 1.8643*x + 1.3886x0.8
4、,1.0时，S =2.5893*x3 - 7.7679*x2 + 6.3643*x - 0.80571输出图如下2-3（1）程序：clear;clc;x1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8;%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1;c=1;for i=1:nc=conv(c,a(i,:);endq=zeros(n,n);r=zeros(n,n+1);for i=1:nq(i,:),r(i,:)=deconv(c,a(i,:);%wn+1/(x-xk)end三次样条插值曲线4次牛顿插值曲线Dw=zeros
5、(1,n);for i=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i);%系数endp=Dw*q;syms x L8;for i=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*);hold ontitle(8次拉格朗日插值);结果如下：8次拉格朗日插值L8 =- 3.2806e-10*x8 + 6.7127e-8*x7 - 5.4292e-6*x6 + 0.00022297*x5 -
6、 0.0049807*x4 + 0.060429*x3 - 0.38141*x2 + 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/x2-y/x;%微分函数Xn=linspace(a
7、,b,n+1);%将定义域分为 n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解for i=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法公式求解for i=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
8、/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp( 改进欧拉法 四阶经典 R-K法);disp(Xn Yn)结果如下：改进欧拉法 四阶经典 R-K法1 10.99887 0.998850.99577 0.99780.99114 0.996940.98532 0.996340.97857 0.996030.97111 0.996060.96311 0.996450.9547 0.997230.94598 0.998410.93705 10.92798 1.0020.91883 1.00440.90964 1.00730.90045 1.01060.89129 1.0
9、1430.88218 1.01840.87315 1.02290.86421 1.02780.85538 1.03310.84665 1.0388（b）程序：clc;clear;a=0;b=1;%定义域H=0.1 0.025 0.01;%步长y0=1/3;%初值f= (x,y) -50*y+50*x2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.2;%准确解Ym=zeros(1,11);for j=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为 n等份Yn=zeros(1,n);%结果矩阵
10、Yn(1)=y0;%赋初值for i=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;endfor k=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
11、四阶经典 R-K法 准确解 误差0.33333 0.33333 04.6055 0.012246 4.593263.062 0.040015 63.022864.05 0.09 863.9611844 0.16 118431.6235e+05 0.25 1.6235e+052.2256e+06 0.36 2.2256e+063.0509e+07 0.49 3.0509e+074.1823e+08 0.64 4.1823e+085.7333e+09 0.81 5.7333e+097.8594e+10 1 7.8594e+10步长为： 2.500000e-02 四阶经典 R-K法 准确解 误差0.
12、33333 0.33333 00.013015 0.012246 0.000768940.040063 0.040015 4.82e-050.090037 0.09 3.6857e-050.16004 0.16 3.6723e-050.25004 0.25 3.6722e-050.36004 0.36 3.6722e-050.49004 0.49 3.6722e-050.64004 0.64 3.6722e-050.81004 0.81 3.6722e-051 1 3.6722e-05步长为： 1.000000e-02 四阶经典 R-K法 准确解 误差0.33333 0.33333 00.01
13、2256 0.012246 9.5673e-060.040016 0.040015 7.8252e-070.090001 0.09 6.6347e-070.16 0.16 6.6226e-070.25 0.25 6.6225e-070.36 0.36 6.6225e-070.49 0.49 6.6225e-070.64 0.64 6.6225e-070.81 0.81 6.6225e-071 1 6.6225e-07由结果可知，步长越小，结果越精确。9-3主程序：clc;clear;options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-4);X,Y
14、= ode45(rigid,0 0.004,1 1 0,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,0 0.004,1 1 0,options);subplot(1,2,2);plot(X,YY(:,1),-,X,YY(:,2),-.,X,YY(:,3),.);title(梯形法);legend(y1, y2,y3);rigid函数function dy = rigid (,y)dy = zeros(3,1); % a column vectordy(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);结果如图：
