东南大学-数值分析上机题作业-MATLAB版

2015.1.9上机作业题报告 1.Chapter 11.1题目设SN=j=2N1j2-1，其精确值为。（1）编制按从大到小的顺序，计算SN的通用程序。（2）编制按从小到大的顺序，计算SN的通用程序。（3）按两种顺序分别计算,并指出有效位数。（编制程序时用单精度）（4）通过本次上机题，你明白了什么？clear;N=input(请输入N值:);Ac=single(3/2-1/N-1/(N+1)/2);Snl2s=single(0);Sns2l=single(0);for i=2:N Snl2s=Snl2s+1/(i*i-1);endfor i=N:-1:2 Sns2l=Sns2l+1/(i*i-1);endfprintf(精确值为： %fn,Ac);fprintf(从大到小的顺序累加得SN=%fn,Snl2s);fprintf(从小到大的顺序累加得SN=%fn,Sns2l);disp(=);1.2程序 P20T17请输入N值:102精确值为： 0.740049从大到小的顺序累加得SN=0.740049从小到大的顺序累加得SN=0.740050= P20T17请输入N值:104精确值为： 0.749900从大到小的顺序累加得SN=0.7498521.3运行结果从小到大的顺序累加得SN=0.749900= P20T17请输入N值:106精确值为： 0.749999从大到小的顺序累加得SN=0.749852从小到大的顺序累加得SN=0.749999=1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。按从小到大的顺序，有效位数分别为：5，6，6。可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。2.Chapter 22.1题目（1）给定初值及容许误差，编制牛顿法解方程f(x)=0的通用程序。（2）给定方程,易知其有三个根由牛顿方法的局部收敛性可知存在当时，Newton迭代序列收敛于根x2*。试确定尽可能大的。试取若干初始值，观察当时Newton序列的收敛性以及收敛于哪一个根。（3）通过本上机题，你明白了什么？2.2程序f(x)函数m文件：fu.mfunction Fu=fu(x)Fu=x3/3-x;endf(x)函数m文件：dfu.mfunction Fu=dfu(x)Fu=x2-1;end用Newton法求根的通用程序Newton.mclear;x0=input(请输入初值x0:);ep=input(请输入容许误差:);flag=1;while flag=1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)epflag=0; endx0=x1;endfprintf(方程的一个近似解为：%fn,x0);寻找最大值的程序：Find.mcleareps=input(请输入搜索精度：);ep=input(请输入容许误差：);flag=1;k=0;x0=0;while flag=1 sigma=k*eps;x0=sigma; k=k+1; m=0; flag1=1; while flag1=1 & m=103 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; endendfprintf(最大的sigma值为：%fn,sigma);2.3运行结果（1）寻找最大的值。算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。运行Find.m，得到在不同的搜索精度下的最大sigma值。 Find请输入搜索精度：10-6请输入容许误差：10-6最大的sigma值为：0.774597 Find请输入搜索精度：10-4请输入容许误差：10-6最大的sigma值为：0.774600 Find请输入搜索精度：10-2请输入容许误差：10-6最大的sigma值为：0.780000（2）运行Newton.m在（-，-1）内取初值，运行结果如下：X0Xk-1000-1.732051-500-1.732051-100-1.732051-10-1.732051-5-1.732051-2.5-1.732051-1.5-1.732051可见，在（-，-1）区间内取初值，Newton序列收敛，且收敛于根-3。在（-1，-）内取初值，运行结果如下：X0Xk-0.951.732051-0.851.732051-0.81.732051-0.7745981.732051可见，在（-1，-）内取初值，Newton序列收敛，且收敛于根3。在（-，）内内取初值，运行结果如下：X0Xk-0.7745960.000000-0.550.000000-0.350.000000-0.150.0000000.050.0000000.250.0000000.450.0000000.650.0000000.7745960.000000可见，在（-，）内取初值，Newton序列收敛，且收敛于根0。在（，1）内取初值，运行结果如下：X0Xk0.774598-1.7320510.8-1.7320510.85-1.7320510.95-1.732051可见，在（，1）内取初值，Newton序列收敛，且收敛于根-3在（1，+）内取初值，运行结果如下：X0Xk1.51.7320512.51.73205151.732051101.7320511001.7320515001.73205110001.732051可见，在（1，+）内取初值，Newton序列收敛，且收敛于根33.Chapter 33.1题目对于某电路的分析，归结为求解线性方程组RI=V，其中（1）编制解n阶线性方程组Ax=b的列主元高斯消去法的通用程序；（2）用所编程序线性方程组RI=V，并打印出解向量，保留5位有效数字；（3）本题编程之中，你提高了哪些编程能力？3.2程序n=input(请输入线性方程组阶数: n=);b=zeros(1,n);A=input(请输入系数矩阵：A=n);b(1,:)=input(请输入线性方程组右端向量：b=n);b=b; C=A,b;for i=1:n-1 maximum,index=max(abs(C(i:n,i); index=index+i-1; T=C(index,:); C(index,:)=C(i,:); C(i,:)=T; for k=i+1:n if C(k,i)=0 C(k,:)=C(k,:)-C(k,i)/C(i,i)*C(i,:); end endend%回代求解x=zeros(n,1);x(n)=C(n,n+1)/C(n,n);for i=n-1:-1:1 x(i)=(C(i,n+1)-C(i,i+1:n)*x(i+1:n,1)/C(i,i);enddisp(方程组的解为:);fprintf(%.5gn,x);3.3运行结果运行程序，输入系数矩阵和方程组右端列向量。运行过程与结果如下图所示： P126T39请输入线性方程组阶数: n=4请输入系数矩阵：A=136.01 90.86 0 0;90.86 98.81 -67.59 0;0 -67.59 132.01 46.26 ;0 0 46.26 177.17请输入线性方程组右端向量：b=-33.254 49.79 28.067 -7.324方程组的解为:-2957.44426.62495-651.49 P126T39请输入线性方程组阶数: n=9请输入系数矩阵：A=31 -13 0 0 0 -10 0 0 0;-13 35 -9 0 -11 0 0 0 0;0 -9 31 -10 0 0 0 0 0;0 0 -10 79 -30 0 0 0 -9;0 0 0 -30 57 -7 0 -5 0;0 0 0 0 -7 47 -30 0 0;0 0 0 0 0 -30 41 0 0;0 0 0 0 -5 0 0 27 -2;0 0 0 -9 0 0 0 -2 29请输入线性方程组右端向量：b=-15 27 -23 0 -20 12 -7 7 10方程组的解为:-0.289230.34544-0.71281-0.22061-0.43040.15431-0.0578230.201050.29023可看出，算得的该线性方程组的解向量为：-0.28923 0.34544 -0.71281 -0.22061 -0.4304 0.15431 -0.057823 0.20105 0.290234.Chapter 44.1题目（1）编制求第一型3次样条插值函数的通用程序；（2）已知汽车门曲线型值点的数据如下：i012345678910Xi012345678910Yi2.513.304.044.705.225.545.785.405.575.705.80端点条件为y0=0.8， y10=0.2，用所编程序求车门的3次样条插值函数S(x),并打印出S(i+0.5)，i=0,1，9。4.2程序cleardigits(6);n=input(请输入节点数：n=);xn=zeros(1,n);yn=zeros(1,n);xn(1,:)=input(请输入节点坐标：);yn(1,:)=input(请输入节点处函数值：);dy0=input(请输入左边界条件：y（x0）=);dyn=input(请输入右边界条件：y（xn）=);%=求d=%d=zeros(n,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n-1 h(i)=xn(i+1)-xn(i); f1(i)=(yn(i+1)-yn(i)/h(i);endfor i=2:n-1 f2(i)=(f1(i)-f1(i-1)/(xn(i+1)-xn(i-1); d(i)=6*f2(i);endd(i)=6*(f1(1)-dy0)/h(1);d(n)=6*(dyn-f1(n-1)/h(n-1);%=求Mi=%A=zeros(n);miu=zeros(1,n-2);lamda=zeros(1,n-2);for i=1:n-2 miu(i)=h(i)/(h(i)+h(i+1); lamda(i)=1-miu(i);endA(1,2)=1;A(n,n-1)=1;for i=1:n A(i,i)=2;endfor i=2:n-1 A(i,i-1)=miu(i-1); A(i,i+1)=lamda(i-1);endM=Ad;%=回代求插值函数=%syms x;for i=1:n-1; Sx(i)=collect(yn(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-xn(i)+M(i)/2*(x-xn(i)2+(M(i+1)-M(i)/(6*h(i)*(x-xn(i)3); Sx(i)=vpa(Sx(i),6);endS=zeros(1,n-1);for i=1:n-1 x=xn(i)+0.5; S(i)=yn(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-xn(i)+M(i)/2*(x-xn(i)2+(M(i+1)-M(i)/(6*h(i)*(x-xn(i)3;end%=打印结果=%disp(S(x)=);for i=1:n-1 format short; fprintf( %s (%dx P195T37请输入节点数：n=11请输入节点坐标：0 1 2 3 4 5 6 7 8 9 10请输入节点处函数值：2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80请输入左边界条件：y（x0）=0.8请输入右边界条件：y（xn）=0.2S(x)= 0.79*x + 0.0158344*x2 - 0.0158344*x3 + 2.51 (0x1)= 0.830013*x - 0.0241785*x2 - 0.00249676*x3 + 2.49666 (1x2)= 0.809832*x - 0.0140879*x2 - 0.00417854*x3 + 2.51012 (2x3)= 0.315407*x2 - 0.178653*x - 0.0407891*x3 + 3.4986 (3x4)= 6.9313*x - 1.46208*x2 + 0.107335*x3 - 5.98133 (4x5)= 4.1762*x2 - 21.2601*x - 0.26855*x3 + 41.0043 (5x6)= 53.8449*x - 8.3413*x2 + 0.426866*x3 - 109.206 (6x7)= 6.27011*x2 - 48.435*x - 0.268915*x3 + 129.447 (7x8)= 14.4854*x - 1.59494*x2 + 0.0587951*x3 - 38.3403 (8x9)= 13.2458*x - 1.45831*x2 + 0.053735*x3 - 34.5615 (9x10)=S(i+0.5) i x(i+0.5) S(i+0.5) 1 0.50000 2.90698 2 1.50000 3.67885 3 2.50000 4.38136 4 3.50000 4.98822 5 4.50000 5.38326 6 5.50000 5.72372 7 6.50000 5.59435 8 7.50000 5.43012 9 8.50000 5.65892 10 9.50000 5.731725.Chapter 55.1题目用Romberg求积法计算积分-1111+100x2dx的近似值，要求误差不超过0.510-7。5.2程序%被积函数m文件：fx.mfunction Fx=fx(x)Fx=1/(1+100*x*x);end%Romberg求积法计算积分的通用程序function Romberg()clear;a=input(请输入积分下限：a=);b=input(请输入积分上限：b=);eps=input(请输入允许精度：eps=);%=计算Tn=%function Tn=T(n)Tn=0;h=(b-a)/n;x=zeros(1,n+1);for k=1:n+1 x(k)=a+(k-1)*h;endfor j=1:n Tn=Tn+h*(fx(x(j)+fx(x(j+1)/2;endend%=计算Sn=%function Sn=S(n)Sn=4/3*T(2*n)-1/3*T(n);end%=计算Cn=%function Cn=C(n)Cn=16/15*S(2*n)-1/15*S(n);end%=计算Rn=%function Rn=R(n)Rn=64/63*C(2*n)-1/63*C(n);end%=计算满足允许精度的Rn，并打印输出=%i=1;flag=1;while flag=1if abs(R(2i)-R(2(i-1)/255 Romberg请输入积分下限：a=-1请输入积分上限：b=1请输入允许精度：eps=0.5*10-7该积分的值为：0.29422555.4结果分析手动化简该定积分并最终求得的值为：0.294225534860747，误差限为：3.48610-8，可见，程序完成了计算要求。6.Chapter 66.1题目常微分方程初值问题数值解（1）编制RK4方法的通用程序；（2）编制AB4方法的通用程序（由RK4提供初值）；（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；（5）对于初值问题y=-x2y2y0=3取步长h=0.1，应用（1）-（4）中的四种方法进行计算，并将计算结果和精确解yx=3/(1+x3)作比较；（6）通过本上机题，你能得到哪些结论？6.2程序%f（x，y）函数m文件：fxy.mfunction FXY=fxy(x,y)FXY=-x*x*y*y;end%精确解y（x）函数m文件：fx.mfunction FX=fx(x)FX=3/(1+x*x*x);end%RK4法通用程序function RK4()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=1:N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i); k2=fxy(x(i)+0.5*h,y(i)+0.5*h*k1); k3=fxy(x(i)+0.5*h,y(i)+0.5*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);end disp(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend%AB4法通用程序function AB4()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=1:N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i); k2=fxy(x(i)+0.5*h,y(i)+0.5*h*k1); k3=fxy(x(i)+0.5*h,y(i)+0.5*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);endfor i=4:N-1 y(i+1)=y(i)+h/24*(55*fxy(x(i),y(i)-59*fxy(x(i-1),y(i-1)+37*fxy(x(i-2),y(i-2)-9*fxy(x(i-3),y(i-3);end disp(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend%AB4-AM4预测校正法通用程序function AB4AM4()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=1:N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i); k2=fxy(x(i)+0.5*h,y(i)+0.5*h*k1); k3=fxy(x(i)+0.5*h,y(i)+0.5*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);endfor i=4:N-1 yp(i+1)=y(i)+h/24*(55*fxy(x(i),y(i)-59*fxy(x(i-1),y(i-1)+37*fxy(x(i-2),y(i-2)-9*fxy(x(i-3),y(i-3); y(i+1)=y(i)+h/24*(9*fxy(x(i+1),yp(i+1)+19*fxy(x(i),y(i)-5*fxy(x(i-1),y(i-1)+fxy(x(i-2),y(i-2);end disp(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend%带改进的AB4-AM4预测校正方法function AB4AM4plus()clear;x(1)=input(请输入初始x值：x0=);y(1)=input(请输入初值条件：y（x0）=);N=input(请输入计算步长：N=);h=input(请输入步长：h=);for i=1:N-1 x(i+1)=x(i)+h; k1=fxy(x(i),y(i); k2=fxy(x(i)+0.5*h,y(i)+0.5*h*k1); k3=fxy(x(i)+0.5*h,y(i)+0.5*h*k2); k4=fxy(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);endfor i=4:N-1 yp(i+1)=y(i)+h/24*(55*fxy(x(i),y(i)-59*fxy(x(i-1),y(i-1)+37*fxy(x(i-2),y(i-2)-9*fxy(x(i-3),y(i-3); yc(i+1)=y(i)+h/24*(9*fxy(x(i+1),yp(i+1)+19*fxy(x(i),y(i)-5*fxy(x(i-1),y(i-1)+fxy(x(i-2),y(i-2); y(i+1)=251/270*yc(i+1)+19/270*yp(i+1);end disp(i xi yi y(xi) y(xi)-yi); disp(-);for i=1:N fprintf(%d %f %f %f %fn,i,x(i),y(i),fx(x(i),fx(x(i)-y(i); disp(-);endend6.3运行结果（1）RK4法 RK4请输入初始x值：x0=0请输入初值条件：y（x0）=3请输入计算步数：N=15请输入步长：h=0.1i xi yi y(xi) y(xi)-yi-1 0.000000 3.000000 3.000000 0.000000-2 0.100000 2.997003 2.997003 0.000000-3 0.200000 2.976190 2.976190 0.000000-4 0.300000 2.921129 2.921130 0.000001-5 0.400000 2.819547 2.819549 0.000002-6 0.500000 2.666663 2.666667 0.000003-7 0.600000 2.467100 2.467105 0.000005-8 0.700000 2.233799 2.233805 0.000006-9 0.800000 1.984123 1.984127 0.000004-10 0.900000 1.735107 1.735107 -0.000000-11 1.000000 1.500006 1.500000 -0.000006-12 1.100000 1.287013 1.287001 -0.000011-13 1.200000 1.099722 1.099707 -0.000015-14 1.300000 0.938397 0.938380 -0.000018-15 1.400000 0.801300 0.801282 -0.000018-（2）AB4法 AB4请输入初始x值：x0=0请输入初值条件：y（x0）=3请输入计算步数：N=15请输入步长：h=0.1i xi yi y(xi) y(xi)-yi-1 0.000000 3.000000 3.000000 0.000000-2 0.100000 2.997003 2.997003 0.000000-3 0.200000 2.976190 2.976190 0.000000-4 0.300000 2.921129 2.921130 0.000001-5 0.400000 2.818389 2.819549 0.001160-6 0.500000 2.664672 2.666667 0.001994-7 0.600000 2.465203 2.467105 0.001903-8 0.700000 2.233079 2.233805 0.000726-9 0.800000 1.984951 1.984127 -0.000824-10 0.900000 1.737043 1.735107 -0.001936-11 1.000000 1.502195 1.500000 -0.002195-12 1.100000 1.288763 1.287001 -0.001762-13 1.200000 1.100724 1.099707 -0.001017-14 1.300000 0.938710 0.938380 -0.000331-15 1.400000 0.801135 0.801282 0.000147-（3）AB4-AM4预测校正法 AB4AM4请输入初始x值：x0=0请输入初值条件：y（x0）=3请输入计算步数：N=15请输入步长：h=0.1i xi yi y(xi) y(xi)-yi-1 0.000000 3.000000 3.000000 0.000000-2 0.100000 2.997003 2.997003 0.000000-3 0.200000 2.976190 2.976190 0.000000-4 0.300000 2.921129 2.921130 0.000001-5 0.400000 2.819678 2.819549 -0.000130-6 0.500000 2.666876 2.666667 -0.000209-7 0.600000 2.467252 2.467105 -0.000147-8 0.700000 2.233731 2.233805 0.000073-9 0.800000 1.983787 1.984127 0.000340-10 0.900000 1.734607 1.735107 0.000500-11 1.000000 1.499516 1.500000 0.000484-12 1.100000 1.286657 1.287001 0.000344-13 1.200000 1.099533 1.099707 0.000174-14 1.300000 0.938343 0.938380 0.000037-15 1.400000 0.801327 0.801282 -0.000045-（4）带改进的AB4-AM4预测校正法 AB4AM4plus请输入初始x值：x0=0请输入初值条件：y（x0）=3请输入计算步数：N=15请输入步长：h=0.1i xi yi y(xi) y(xi)-yi-1 0.000000 3.000000 3.000000 0.000000-