工程数值计算matlab试验报告-第三次试验

《工程数值计算》上机实验报告(第三次)学生姓名**** 班级**** 学号*****任课教师*****上机时间2019年n月7日,报告完成2019年n月8日1.实验目的：简述实验任务和目的任务3.1：采用以下两种方法求解线性方程组，2%]+24+3々-34x+7xn+7xo=1TOC\o"1-5"\h\z12 3{—2x+ =-712 3(1)逆矩阵法(2)数值迭代法任务3.2：非线性方程(组)的求根(1)求2炉—4%4+3%3—8x2+9=。的根(2))伞兵降落速度方程为v(t)=即(1一e-U)，其中g=9.81m/s"2,m=68kg，若要下落Ct=10s时速度达到v=40m/s请问阻力系数c是多少?(3)求方程的解x2+xy+y=3%2—4%+3=0任务3.3：微分方程组的求解(3)求方程的解—二(j(y—x)dt已知大气动力学模型的洛伦兹方程^=rx-y-xz，其中x(t)为大气流动的强度,y(t)dt=xy—bzdt和Z(t)为水平和垂直方向的温度处化。如果。=10,b=2.666667,r=28,初始条件为x(0)=6,y(0)=z(0)=5,请你：(1)计算t=0〜20范围内x(t)、y(t)和z(t)的数值解。(2)分别画出y(t)和z(t)随时间t的变化曲线y-t和z-t。(3)分别画出x随y和x随z变化的曲线(相平面图)：x-y和x-z。(4)轻微改变初始条件x(0)=6.01,y(0)=z(0)=5,重复以上过程，画出相平面图(x-y和x-z),并与以上结果进行比较。2.计算方法：针对实验任务，结合课堂内容，说明解决方法，如：采用何种1/10

理论，列出相关公式，说明计算步骤，写出程序框图等任务3.1两种方法求解线性方程组(1)逆矩阵法计算原理:对于a11a21a12a22a1na21X1X2b1b2或对于a11a21a12a22a1na21X1X2b1b2或Ax=bn1n2ann解为a11121na21a22a21b1b2或*=A1bama.n2nn流程图:流程图:(2)数值迭代法(ax+—+ax(2)数值迭代法(ax+—+ax111 1nnQX+^+UX211 .2nnan1X1+-+annXnb2将对角线上各项留在左边，其余移到右边，方程写为1x=——(ax+ax+—+ax+axb)1a 122 133 1n1n11nn1111x=—(ax+ax+ +ax+axb)2ay211 233 2n1n12nn2221,x= (ax+ax+ax+ +ax b)naKn11n22n33 nn1n1 nnn写成迭代形式2/10”1)

1%但)

2-1-1a”22(a0134k)”1)

1%但)

2-1-1a”22(a0134k)(a21X1k1呼),、一1x(k1)= (ax(k1)ax(k1)ax(k1)1aam1 n22 a33nnax(k1)—b)

nnn 1经过多次迭代计算逐渐逼近真值。流程图:任务3.2：非线性方程(组)的求根(1)计算原理：roots(P)函数将P视为一个具有n+1个元素的向量，代表nXn矩阵A的n次特征多项式。多项式的根通过计算伴随矩阵A的特征值得出。已知%%九…a1x%=0,P=[4%一]…%]3/10

a儿一1—a n~2aann1a儿一1—a n~2aann10构造矩阵般"01[0… 一^rran… 0… 0•,••.0 1]设Ax=^x（x为该方程的解）以3项来说明原式：a3x3+a2x2+4%1+a%A=[a310ai-101a0一力00a—2*xa3X川」第3aa 1*% 0*%=Ax

1a2a3 13 3X1=A%,%=AX3联立以上3个方程得:a/3+a/+q/1+%=0因此，只要解出入（特征根）的值，就得到a/3+a/2+a/1+%=0的解。从而，我们把解方程的根转化成了特征根的求解。流程图：（2）定义自定义函数，调用matlab的fslove函数求解方程。（3）调用matlab的slove函数求解方程。任务3.3：微分方程组的求解计算原理：欧拉公式2+1=〃+f（xn,ynM （初值条件｛；；=；=；））流程图：4/10

3.程序设计：根据前面提到的计算方法编写程序；写出程序代码，并结合计算方法对程序中关键步骤进行必要的文字说明任务3.1两种方法求解线性方程组%逆矩阵法求解方程组clear;A=[223;477;-245];%（系数矩阵）B=[3;1;-7];%（列向量）C=inv（A）*B%（求解）%数值迭代法求解n=length（B）;N=20;eb=1e-3;%（迭代收敛精度）5/10%给定迭代初始点fori=1:nx(1,i)=0;endfori=1:Nix(i+1,1)=(B(1)-A(1,2)*x(i,2)-A(1,3)*x(i,3))/A(1,1);%(迭代格式计算)x(i+1,2)=(B(2)-A(2,1)*x(i+1,1)-A(2,3)*x(i,3))/A(2,2);%(迭代格式计算)x(i+1,3)=(B(3)-A(3,1)*x(i+1,1)-A(3,2)*x(i+1,2))/A(3,3);%(迭代格式计算)e(i)=abs(x(i+1,1)-x(i,1));ife(i)<eb%(判断精度)break%(跳出循环)endendx任务3.2：非线性方程(组)的求根%非线性方程组求根(1)clear;c=[2,-4,3,-8,9];%(方程系数)x=roots(c)%定义函数functionF=mf(x)globalmgtvF=m*g/x*(1-exp(-x/m*t))-v;%调用函数globalmgtvm=68;6/10g=9.8;t=10;v=40;X0=12;%(定义初值)[X,Fval,Exit]=fsolve(@mf,12)clear;[x,y]=solve(，x"2+x*y+y-3=0，,，x"2-4*x+3=0，)任务3.3：微分方程组的求解%微分方程求根clear;T=20;%(计算时长)dt=2.5;%(时间步长)N1=T/dt;%初始点定义x=6;y=5;z=5;%参数取值a=10;b=2.666667;c=28;%Part2-程序主体x1(1)=x;y1(1)=y;z1(1)=z;t(1)=0;fori=1:N1t(i+1)=i*dt;7/10x1(i+1)=x1(i)-a*(x1(i)-y1(i))*dt;y1(i+1)=y1(i)+(c*x1(i)-y1(i)-x1(i)*z1(i))*dt;z1(i+1)=z1(i)+(x1(i)*y1(i)-b*z1(i))*dt;end%Part3-画图figure;subplot(2,1,1)plot(t,y1,'r-')subplot(2,1,2)plot(t,z1,'b-')figure;subplot(1,2,1)plot(y1,x1,'r-')subplot(1,2,2)plot(z1,x1,，b-，)4.结果分析：给出计算结果(可用数值、图表、曲线表示)，并进行分析任务3.1两种方法求解线性方程组(1)逆矩阵法2.0000-2.00001.0000(2)数值迭代法8/10

16经过16次迭代计算，得出符合精度需要的值（末尾一行）。经过16次迭代计算，得出符合精度需要的值（末尾一行）。任务3.2：非线性方程（组）的求根001.5000-0.7143-0.22S62.S57L-L.O0Q80.49471.8478-L.40770.46S32.2098-L.S85L0.75201.Q571-L.727S0.76492.0802-L.31070.S8071.Q897-L.87480.SQ572.0312-L.31350.Q132L998694250,9534工012496050,9733[ooo.--L,97360,97912.0050-L.99190,9875工000698790,99062.0020<.99170,99422.0004-L99450,99572.0009-L.99620,9973（1）-0.4675+L.356Gi-0.4675-1.356614675+0..78314675-Q.17S31（2）14.7585Fval=2891e-O8Exit千Exit=1，说明计算结果可信。9/9/10(3)-1-3/2任务3.3：微分方程组的求解原题数据输出结果（左y-t；右z-t）工，而皿 1在———.——,—— |产更改数据后（第四问要求）结果（左x-y；右x-z）更改数据