数值分析实验题答案

数值分析实验报告姓名：院系：能源学院热能工程学号：2014年4月2/8习题一实验3.2编制正交化多项式最小二乘拟合程序，并用于求解3次多项式最小二乘拟合问题，作拟合曲线的图形，计算平方误差，并与以函数0nkkx为基的多项式最小二乘拟合的结果作比较表1ix-1.0-0.50.00.51.01.52.0iy-4.447-0.4520.5510.0480.4470.5494.552首先使用以函数0nkkx为基的多项式最小二乘拟合，代码如下：x=linspace(-1,2,7);y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];fori=1:4forj=1:4A(i,j)=sum(x.^(3+i-j));endB(i)=sum(y.*x.^(i-1));endal=A\B';deltal=sum((polyval(al,x)-y).^2);然后使用正交化多项式方法作最小二乘拟合并画图，代码如下：[P,phy]=orthpoly(x,y,3);delta2=sum((polyval(P,x)-y).^2);plot(x1,polyval(a1,x1),’r-’,x1,polyval(P,x1),’b-’,x,y,’k.’);%orthpolyfuctionfunction[P,phy]=orthpoly(x,y,k)symsX;phy=zeros(k+1,length(x));phy(1,:)=1;P(1)=sum(y)/length(x);a=sum(x.*phy(1,:).*phy(1,:))/sum(phy(1,:).*phy(1,:));phy(2,:)=x-a;Phy(1)=sym(1);Phy(2)=X-a;P(2)=sum(y.*phy(2,:))/sum(phy(2,:).*phy(2,:));fori=3:k+13/8a=sum(x.*phy(i-1,:).*phy(i-1,:))/sum(phy(i-1,:).*phy(i-1,:));b=sum(phy(i-1,:).*phy(i-1,:))/sum(phy(i-2,:).*phy(i-2,:));phy(i,:)=(x-a).*phy(i-1,:)-b.*phy(i-2,:);Phy(i)=(X-a)*Phy(i-1)-b*Phy(i-2);P(i)=sum(y.*phy(i,:))/sum(phy(i,:).*phy(i,:));endP=sym2poly(sum(P.*Phy));end拟合得到的图形如下（图1）：从图形来看，二者与数据点都很吻合。计算结果为：delta1=2.1762e-05delta2=4.4701e-04可以看出对于3次多项式以0nkkx为基的多项式拟合精度更高。图14/8习题二实验4.2分别用复化Simpson公式与变步长Simpson公式计算，要求绝对误差限为71=102，输出每种方法所需的节点数和积分近似值并分析比较结果。（1）6220()10xxxdx（2）10xxdx（3）6220()10xxxdx对于复化Simpson公式，使用事前误差估计法得到所需计算节点数，有以下误差估计公式：4(4)()()()(),(,)1802nbahRffab对（1）式有：(4)22max()36144xfxx0.0266h75.2baNh取76N，计算节点数为21153NN对（2）（3）式由于其4阶导数值分布极不均匀，用最大值来估计所需计算节点数造成很大浪费，尝试多次后分别取153N和1091N代码如下：formatlongf=@(x)x.^6/10-x.^2+x;disp(’4.2计算结果_1’);Simpson(f,0,2,153);fprintf(’标准值=%.9f\n\n’,128/70-8/3+2);f=@(x)x.^1.5;disp(’4.2计算结果_2’);Simpson(f,0,1,153);fprintf(’标准值=%.9f\n\n’,0.4);f=@(x)x.^(-0.5);disp(’4.2计算结果_3’);Simpson(f,5,200,1091);fprintf(’标准值=%.9f\n\n’,2*(sqrt(200)-sqrt(5)));%SimpsonfunctionFunctionSimpson(f,a,b,N)s=@(f,x,N)(x(N)-x(1))/(N-1)*2/6*(f(x(1))+f(x(N))+4*sum(f(x(2:2:N-1)))+2*sum(f(x(3:2:N-2))));x=linspace(a,b,N);s1=s(f,x,N);[t1,N1]=quad(f,a,b,0.5e-7);disp(’用复化公式计算:Simpson’);fprintf(’节点数:%d/n近似值:%.9f\n’,N,s1);disp(’用变步长公式计算:Simpson’)5/8fprintf(’节点数:%d/n近似值:%.9f\n’,N1,t1);end计算结果如下：（1）计算结果用复化Simpson公式计算：节点数：153近似值：1.161904777用变步长Simpson公式计算：节点数：77近似值：1.161904766标准值：1.161904762（2）计算结果用复化Simpson公式计算：节点数：153近似值：0.400000049用变步长Simpson公式计算：节点数：33近似值：0.400000069标准值：0.400000000（3）计算结果用复化Simpson公式计算：节点数：1091近似值：23.812135331用变步长Simpson公式计算：节点数：145近似值：23.812135297标准值：23.812135292从以上计算结果可以看出，变步长simpson公式所需节点数明显减少，因为3个函数的4阶导数在积分区间内分布都部均匀，（2）（3）式更为严重，在先范围区间内导数值远大于其他区间，只需要在这些区间增加节点数就可以达到指定精度，而Simpson公式需要在全体积分限内采用较小间距才满足条件。6/8习题三实验5.3常微分方程初值问题{𝑦′=−𝑦+𝑐𝑜𝑠2𝑥−2𝑠𝑖𝑛2𝑥+2𝑥𝑒−𝑥0𝑥2𝑦(0)=1有精确解2()cos2xyxxex(1)选择一个步长h使4阶Adams预测校正修正法和经典R-K法均稳定，分别用这两种方法计算初值问题，以表格形式列出10个等距节点上的计算值和精确值并比较他们的计算精度。(2)取h=0.001，依然用这两种方法计算，比较两种方法所花的计算机CPU时间。代码如下：rhs=@(x,y)-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x);std=@(x)x.*x.*exp(-x)+cos(2*x);h=0.2;x=0:h:2;y1=rk4(rhs,0,1,h,2);y2=adams_pc(rhs,0,1,h,2);disp([x;y1;y2;std(x)]);plot(x,abs((y1-std(x))./std(x)),’r.’,x,abs((y2-std(x))./std(x)),’b*’)h=0.001;x=0:h:2;tic;y1=rk4(rhs,0,1,h,2);time1=toc;tic;y2=adams_pc(rhs,0,1,h,2);time2=toc;disp([x(length(x));y1(length(x));y2(length(x));std(x(length(x)))]);%rk4functionfunctiony=rk4(rhs,x0,y0,h,x1)x=x0:h:x1;y=zeros(size(x));y(1)=y0;fori=1:length(x)-1k1=rhs(x(i),y(i));k2=rhs(x(i)+h/2,y(i)+h/2*k1);k3=rhs(x(i)+h/2,y(i)+h/2*k2);k4=rhs(x(i)+h,y(i)+h*k3);y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);endend7/8%Adamsfunctionfunctiony=adams_pc(rhs,x0,y0,h,x1)x=x0:h:x1;y=zeros(size(x));rhs_tmp=zeros(size(x));y=rk4(rhs,x0,y0,h,3*h);ync=0;ynp=0;rhs_tmp(1)=rhs(x(1),y(1));rhs_tmp(2)=rhs(x(2),y(2));rhs_tmp(3)=rhs(x(3),y(3));fori=4:length(x)-1rhs_tmp(i)=rhs(x(i),y(i));yp=y(i)+h/24*(55*rhs_tmp(i)-59*rhs_tmp(i-1)+37*rhs_tmp(i-2)-9*rhs_tmp(i-3));ypm=yp+251/270*(ync-ynp);ync=y(i)+h/24*(9*rhs(x(i+1),ypm)+19*rhs_tmp(i)-5*rhs_tmp(i-1)+rhs_tmp(i-2));y(i+1)=ync-19/270*(ync-yp);ynp=yp;endend输出的结果如下：节点x坐标经典R-K方法Adams预估-校正修正法标准值0.01.0000000000000001.0000000000000001.0000000000000000.20.9538097986034620.9538097986034620.9538102241260040.40.8039571873024350.8039571873024350.8039579167128680.60.5599301449143590.5599301449143590.5599299434705230.80.2583736455298820.2580864435878330.2583710147337331.0-0.048261190980325-0.048506026987456-0.0482673953757001.2-0.303663865842853-0.303846236688374-0.3036740503876741.4-0.458878611459712-0.458961883351702-0.4588922913431091.6-0.481423834832514-0.481371870814424-0.4814396897284351.8-0.361173917672502-0.360982088404667-0.3611900184962072.0-0.112288336571588-0.111970963792309-0.112302487917161对于h=0.2，Adams方法的误差要高于经典R-K方法。h=0.001时，运行结果为：节点x坐标经典R-K方法Adams预估-校正修正法标准值2.0-0.112302487917154-0.112302487917161-0.112302487917161经典R-K方法运行时间为0.008448s，Adams方法运行时间为0.009849s，速8/8度基本持平，Adams方法在积分终点的值更精确。若继续减小h=0.0001，两者运行时间分别为0.077s和0.254s，Adams方法速度明显慢于经典R-K方法。