数值分析课程设计报告.doc

数值分析课程设计报告 院 系 理学院 班 级 11级应用数学 学生姓名 汪佩祺 学 号 201130512276 任课教师 雷秀仁 完成日期 2013年8月23日 目录设计题一（1） 问题 1（2） 设计思路 1（3） 程序清单 1（4） 程序运行与结果 5（5） 结果分析 6设计题二（1）问题 7（2）设计思路 7（3）程序清单 7（4）程序运行与结果 9（5）结果分析 11设计题四（1）问题 11（2）设计思路 12（3）程序清单 12（4）程序运行与结果 12（5）结果分析 13设计题五（1）问题 14（2）设计思路 14（3）程序清单 14（4）程序运行与结果 20（5）结果分析 24心得体会 25自我评价 26设计题一1） 问题编写解线性代数方程组的列主元高斯消去法的函数，并调用该函数计算某个9阶以上的非奇异阵A的逆矩阵。通过计算AA-1检查答案，并与使用inv（A）所得结果和运行时间进行比较。2） 设计思路分为选主元、消去、回代和矩阵求逆四个过程：对于第一列，先选择最大的元素作为主元，并将该行于第一行对换，再进行消去过程。对于第k列，对比由k到n的元素最大值作为主元，并将该行与第k行对换，再进行消去过程。最终，该行列式转化为一个上三角阵，紧接着进行回代过程。3） 程序清单M文件：%Gaosixiaoqu.mfunctionx=Gaosixiaoqu(a,b)n=length(a);x=zeros(n,1); a=a b;for m=1:n-1 max=m; for i=m+1:n %选取主元 if a(i,m)a(max,m) max=i; end end temp=a(m,m:n+1); a(m,m:n+1)=a(max,m:n+1); a(max,m:n+1)=temp; for i=m+1:n %消去过程 a(i,m)=-a(i,m)/a(m,m); a(i,m+1:n+1)=a(i,m+1:n+1)+a(i,m)*a(m,m+1:n+1); end end x(n,1)=a(n,n+1)/a(n,n); %回代过程 for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+x(j,1)*a(i,j); end x(i,1)=(a(i,n+1)-sum)/a(i,i); end命令文件： a=1 3 2 9 1 4 8 7 6;2 3 4 9 7 0 1 8 1;2 7 5 0 9 8 9 2 4;3 2 1 0 4 4 6 2 7;1 2 0 9 3 8 5 4 1;2 2 3 4 5 5 1 2 0;5 1 2 3 4 5 6 7 3;1 1 9 2 3 4 5 1 2;2 3 4 1 2 7 4 1 3 b=2;5;0;3;1;2;3;4;8 Gaosixiaoqu(a,b)ans = -14.6399 -5.9055 1.2273 -5.9630 2.4398 7.2287 -6.3115 10.61505.0669 b1=1;0;0;0;0;0;0;0;0 x1=Gaosixiaoqu(a,b1)x1 = 2.1040 1.0097 -0.0492 0.8453 -0.5371 -0.9420 0.6956 -1.3804 -0.5398 b2=0;1;0;0;0;0;0;0;0 x2=Gaosixiaoqu(a,b2)x2 = -1.7759 -0.7142 0.0719 -0.6657 0.4042 0.7226 -0.6298 1.2118 0.5046.x=x1:x2:x3:x4:x5:x6:x7:x8:x9ans = 2.1040 -1.7759 0.1792 -0.6625 -1.9940 4.1295 0.3231 -0.3950 -1.7044 1.0097 -0.7142 0.1930 -0.4288 -0.9521 1.7366 0.0459 -0.3053 -0.5631 -0.0492 0.0719 -0.0321 -0.0187 -0.0043 -0.1076 -0.0154 0.0999 0.1110 0.8453 -0.6657 0.0458 -0.2217 -0.6855 1.6136 0.0161 -0.1118 -0.7253 -0.5371 0.4042 -0.0124 0.2541 0.4710 -0.8117 -0.0770 0.1139 0.2073 -0.9420 0.7226 -0.1287 0.2818 0.8782 -1.6985 -0.0860 0.1296 0.8641 0.6956 -0.6298 0.1619 -0.2518 -0.5503 1.1918 0.1342 -0.0366 -0.7359 -1.3804 1.2118 -0.1440 0.3492 1.2259 -2.7525 -0.0378 0.1838 1.2409 -0.5398 0.5046 -0.1444 0.3619 0.4777 -1.0513 -0.1605 0.0820 0.5395 inv(a)ans = 2.1040 -1.7759 0.1792 -0.6625 -1.9940 4.1295 0.3231 -0.3950 -1.7044 1.0097 -0.7142 0.1930 -0.4288 -0.9521 1.7366 0.0459 -0.3053 -0.5631 -0.0492 0.0719 -0.0321 -0.0187 -0.0043 -0.1076 -0.0154 0.0999 0.1110 0.8453 -0.6657 0.0458 -0.2217 -0.6855 1.6136 0.0161 -0.1118 -0.7253 -0.5371 0.4042 -0.0124 0.2541 0.4710 -0.8117 -0.0770 0.1139 0.2073 -0.9420 0.7226 -0.1287 0.2818 0.8782 -1.6985 -0.0860 0.1296 0.8641 0.6956 -0.6298 0.1619 -0.2518 -0.5503 1.1918 0.1342 -0.0366 -0.7359 -1.3804 1.2118 -0.1440 0.3492 1.2259 -2.7525 -0.0378 0.1838 1.2409 -0.5398 0.5046 -0.1444 0.3619 0.4777 -1.0513 -0.1605 0.0820 0.5395验证得x=inv(a)，二者结果一致。4） 程序运行与结果 M文件过程图1 过程图2结果图5）结果分析本程序利用列主元消去法解矩阵的逆主要是通过解AX=b，根据矩阵求逆过程中的初等变换为依据，将单位矩阵分解成n个列向量，依次求出AXi=bi的解，从而得到inv（a）=x1:x2:x3:x4:x5:x6:x7:x8:x9，经过验证与结果完全一致。cond(a)= 62.4754，是个大数，故a，b的扰动引起的解x的相对误差就比较大，所以稳定性比较弱。优点这个程序过程清晰易懂，明确的分清了换行、消去过程、回代过程等，步骤明确，将高斯法解线性方程组运用到矩阵的求逆上，是对书本知识的一个应用之一。缺点但是，求逆的过程中没有使用较好的方法，浪费了时间。设计题二1） 问题 对于迭代法， 它显然有不动点。 试不用判定收敛阶的定理，设计1至2个数值实验（其中必须有一个不是直接用收敛阶的定义）得到收敛阶数的大概数值。2） 设计思路I、直接用定义求收敛阶：迭代过程显然有不动点。记迭代误差，如果存在实数p=1和非零常数c，使得，则迭代过程是p阶收敛的。II、间接用定义求收敛阶：通过分析题中迭代法，可以知道只需判断ek+1的值即可，因为ek在非收敛阶情况下，去极限值恒为0。3） 程序清单m文件I、直接用定义求收敛阶：%Qiujie.mfunction n=Qiujie(x0)m=0;syms x;x1=0.99*x-x2;k=(abs(x0-x1)/(abs(x0-x)m;while limit(k,x,x0)=0 m=m+1; k=(abs(x0-x1)/(abs(x0-x)m; x=x1;endn=m;II、间接用定义求收敛阶：%Aqiujie.mfunction =Aqiujie()x=0; e=x2-0.99*x;ek=-x;m=0;while abs(e)=0 %判断条件 m=m+1; e=diff(e); enddisp(收敛阶数为p=); %输出收敛阶和收敛的性质disp(p);if m=1 disp(迭代法线性收敛);elseif m=1 disp(迭代法超线性收敛)endend命令文件方法一： Qiujie(0)ans = 1方法二： Aqiujie收敛阶数为p= 1迭代法线性收敛4）程序运行与结果 方法一m文件方法二m文件方法一运行结果方法二运行结果5） 结果分析方法一： 该收敛公式对于x=0是1阶收敛的。方法二：从输出结果可以知道，这个迭代法在不动点x=0处是1阶收敛的。优点一、提供了两种满足题目要求的方法。二、程序输出都有相对详细的提示语，可以比较清楚地看到输出的是几阶收敛。缺点方案具有明显的局限性，仅仅是围绕定义来设计。 设计题四1）问题某飞机头部的光滑外形曲线的型值点坐标由下表给出:012345678910070130210337578776101211421462184105778103135182214244256272275试建立其合适的模拟曲线（未必是用拟合方法），并求在点x100，250，400，500，800处的函数值y及一阶、二阶导数值y，y”。绘出模拟曲线的图形。2）设计思路先根据给定的数据用画图指令画出大致图像，分析其大概形状，然后建立一条比较合适的曲线去逼近原函数。3） 程序清单format longx1=0 70 130 210 337 578 776 1012 1142 1462 1841;y1=0 57 78 103 135 182 214 244 256 272 275;x0=100 250 400 500 800;plot(x1,y1)polyfit(x1,y1,2); %拟合出给出数据函数y2=-0.0001*x1.2+0.3245*x1+28.0193;nihehanshu(x0)plot(x1,y1,b,x1,y2,r)title(模拟曲线和数据曲线)xlabel(x); ylabel(y);legend(数据曲线,拟合曲线);grid on;4）程序运行与结果点x100，250，400，500，800处的函数值y及一阶、二阶导数值y，y”如下：x100250400500800y59.4693102.8943141.8193165.2693223.6193y0.1138-0.2019-0.5178-0.7284-1.3602y”-0.002-0.0021-0.0021-0.0021-0.0025） 结果分析通过Matlab可算得一元二次函数的的表达式：y=-0.0001*x2+0.3245*x+28.0193。优点 I、程序简单易懂。II、经过对比，拟合结果曲线有比较好的拟合，结果有一定程度的可靠性。缺点 程序过于简单，没有深入思考。结果不一定是最佳的拟合曲线。 设计题五1） 问题给定初值问题 其精确解为 ，分别按下列方案求它在节点 处的数值解及误差。比较各方法的优缺点，并将计算结果与精确解做比较（列表、画图）。（方案I）欧拉法，步长h = 0.025, h = 0.1；（方案II）改进的欧拉法，步长h = 0.05, h = 0.1；（方案III）四阶经典龙格库塔法，步长h = 0.1。2） 设计思路1、 设计出题目要求的三种方案的m文件以及微分方程的m文件。2、 代入初值得到各方法的数值解。3、 对比精确解，进行误差分析。4、 作图进行对比。3） 程序清单微分方程m文件%fun.mfunction z=fun(x,y)z=2/x*y+x.2*exp(x);欧拉法m文件%Euler.m function x=Euler(f,x0,y0,xn,M)%f为一阶常微分方程的一般表达式的右端函数;x0，y0 为初始条件;xn 为取值范围的一个端点;M为区间的个数x=zeros(1,M+1); %x为Xn构成的向量y=zeros(1,M+1); %y为Yn构成的向量x(1)=x0;y(1)=y0;h=(xn-x0)/M;for n=1:M x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(f,x(n),y(n);endT=x,y改进欧拉法m文件%GaijinEuler.m function x=GaijinEuler(f,x0,y0,xn,M)x=zeros(1,M+1); %x为Xn构成的向量y=zeros(1,M+1); %y为Yn构成的向量x(1)=x0;y(1)=y0;h=(xn-x0)/M;for n=1:M; x(n+1)=x(n)+h; z0=y(n)+h*feval(f,x(n),y(n); y(n+1)=y(n)+h/2*(feval(f,x(n),y(n)+feval(f,x(n+1),z0);endT=x,y四阶经典龙格库塔法m文件%FRonK.mfunction R=FRonK(f,x0,y0,xn,M);h=(y0-x0)/M;x=zeros(1,M+1); y=zeros(1,M+1); T=x0:h:y0;y(1)=xn;for n=1:M k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+h/2,y(n)+k1/2); k3=h*feval(f,T(n)+h/2,y(n)+k2/2); k4=h*feval(f,T(n)+h,y(n)+k3); y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6;endR=T,y命令文件欧拉法 h=0.025时 clear all; x0=1;y0=0;xn=2;M=40; x=Euler(fun,x0,y0,xn,M) T = 1.0000 0 1.0250 0.0680 1.0500 0.1445 1.0750 0.2301 1.1000 0.3255 1.1250 0.4311 1.1500 0.5478 1.1750 0.6760 1.2000 0.8165 1.2250 0.9701 1.2500 1.1374 1.2750 1.3192 1.3000 1.5164 1.3250 1.7297 1.3500 1.9601 1.3750 2.2085 1.4000 2.4757 1.4250 2.7629 1.4500 3.0709 1.4750 3.4009 1.5000 3.7539 1.5250 4.1311 1.5500 4.5337 1.5750 4.9630 1.6000 5.4201 1.6250 5.9065 1.6500 6.4235 1.6750 6.9725 1.7000 7.5551 1.7250 8.1728 1.7500 8.8272 1.7750 9.5200 1.8000 10.2529 1.8250 11.0277 1.8500 11.8464 1.8750 12.7107 1.9000 13.6228 1.9250 14.5847 1.9500 15.5985 1.9750 16.66672.0000 17.7914x = Columns 1 through 14 1.0000 1.0250 1.0500 1.0750 1.1000 1.1250 1.1500 1.1750 1.2000 1.2250 1.2500 1.2750 1.3000 1.3250 Columns 15 through 28 1.3500 1.3750 1.4000 1.4250 1.4500 1.4750 1.5000 1.5250 1.5500 1.5750 1.6000 1.6250 1.6500 1.6750 Columns 29 through 41 1.7000 1.7250 1.7500 1.7750 1.8000 1.8250 1.8500 1.8750 1.9000 1.9250 1.9500 1.9750 2.0000h=0.1时 clear all; x0=1;y0=0;xn=2;M=10; x=Euler(fun,x0,y0,xn,M) T = 1.0000 0 1.1000 0.2718 1.2000 0.6848 1.3000 1.2770 1.4000 2.0935 1.5000 3.1874 1.6000 4.6208 1.7000 6.4664 1.8000 8.8091 1.9000 11.74802.0000 15.3982x = 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000改进欧拉法 h=0.05时 clear all; x0=1;y0=0;xn=2;M=20; x=GaijinEuler(fun,x0,y0,xn,M) T = 1.0000 0 1.0500 0.1532 1.1000 0.3449 1.1500 0.5801 1.2000 0.8643 1.2500 1.2032 1.3000 1.6031 1.3500 2.0710 1.4000 2.6142 1.4500 3.2406 1.5000 3.9589 1.5500 4.7785 1.6000 5.7092 1.6500 6.7621 1.7000 7.9487 1.7500 9.2816 1.8000 10.7744 1.8500 12.4418 1.9000 14.2993 1.9500 16.3641 2.0000 18.6542 h=0.1时 clear all; x0=1;y0=0;xn=2;M=10; x=GaijinEuler(fun,x0,y0,xn,M） T = 1.0000 0 1.1000 0.3424 1.2000 0.8583 1.3000 1.5927 1.4000 2.5983 1.5000 3.9364 1.6000 5.6789 1.7000 7.9092 1.8000 10.7245 1.9000 14.2374 2.0000 18.5789四阶经典龙格库塔法 h=0.1时 clear all; x0=1;y0=0;xn=2;M=10; R=FRonK(fun,x0,y0,xn,M) R = 1.0000 2.0000 0.9000 1.4105 0.8000 0.9647 0.7000 0.6348 0.6000 0.3975 0.5000 0.2327 0.4000 0.1239 0.3000 0.0570 0.2000 0.0205 0.1000 0.0049 0 Inf画图程序 精确解函数m文件 %f.m function z=f(x);z=x.2*(exp(x)-exp(1); 画图命令 x=1.1:0.1:2x = 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y=f(1.1) f(1.2) f(1.3) f(1.4) f(1.5) f(1.6) f(1.7) f(1.8) f(1.9) f(2)y = 0.3459 0.8666 1.6072 2.6204 3.9677 5.7210 7.9639 10.7936 14.3231 18.6831 plot(x,y).4）程序运行与结果 欧拉法m文件改进欧拉法m文件四阶经典龙格库塔法m文件 部分运行截图 精确解函数图像 欧拉法图像（h=0.1）改进