数值分析大作业2014

1、eg第流Z久大枣课程设计课程名称：高等数值计算设计题目：数值计算B课程设计学号：姓名：完成时间：2014年10月20日题目一：非线性方程求根用Newton法计算下列方程3(1) XX1=0,初值分别为x0=1,%=0.45,Xo=0.65；32(2) x+94x-389X+294=0其三个根分别为13-98。当选择初值xo=2时给6出结果并分析现象，当名=5父10,迭代停止。一、摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。本实验通过使用常用的求解方法二分法和Newton法及改进的Newton法处理几个题目，分析并总结不同方法处理问题的优缺点。观察迭代次数，收敛速度及初

2、值选取对迭代的影响。二、数学原理构造迭代函数的一条很重要的途径是，用近似方程来代替原方程去求根。因此，如果能将非线性方程用线性方程来代替的话，求近似根问题就很容易解决，而且十分方便。Newton法就是把非线性方程线性化的一种方法。在求解非线性方程f(x)=o时，它的困难在于f(x)是非线性函数，为克服这一困难，考虑它的线性展开。设当前点为xk,在xk处的Taylor展开式为f(x)：.f(xk)f(xk)(x-xk)令f(x)=0,可以得到上式的近似方程f(xk)-f(xk)(x-xk)=0设f(xk)#0,解其方程得到xk1=xk-f(xk)(k=。1)f(xk)这就是牛顿迭代公式。用牛顿迭

3、代公式求方程f(x)=0根的方法称为牛顿迭代法。牛顿迭代法的几何意义为，不断用切线来近似曲线得到方程的根，我们知道方程f(x)=0的实根x*是函数y=f(x)的图形与横坐标的交点，xk+是函数f(x)在点(xk,f(xk)处的切线与x轴的交点，此时就是用切线的零点代替曲线的零点，因此，牛顿迭代法又称为切线法。三、程序设计基于MATLAB软件编写程序，先定义一个用Newton法求解的功能函数，然后调用函数用于计算不同的方程。各变量定义见程序。1、选取初值。2、利用公式求解X-=Xk_x4(k=0,1,)f(Xk)3、计算所得Xk+-Xk是否满足精度要求4、如不满足继续迭代运算，如满足则输出所求结

4、果四、结果分析和讨论1、第一题计算结果：首先得到函数在区间卜2.5,2.5的图像，即可知函数f(X)与x轴有交点，也就是说有根，并且从图中能够大致估算到根的位置C口加亦&iidWindow(1)、取初值X。=1时得到卞g值r=1.3247,迭代次数t=4次(2)、取初值X。=0.45时得到卞g值r=1.3247,迭代次数t=42次(3)、取初值X=0.65时得到卞g值r=1.3247,迭代次数t=8次根据结果可以分析得到，当使用牛顿迭代法时，所选初始值对迭代速度(迭代次数)有较大影响。当初始值X。充分接近方程的单根时，可保证迭代序列快速收敛，当初值选择不当时会造成迭代次数大幅增加或不一

5、定收敛。2、第二题计算结果：初值X。=2时，得到根值r=-98,迭代次数为1次。>>工整M.t=HevtanRoct?('x'3-!f-1fc)root=L师»root,tj=KertonfcMrtS(-x-3-z-l'Fth4S)工力L受打t42»rwtatl=Kswt口由ootZt'k3rT：.必651rocrt=L32dT1=A3»cootjt=JIewtonRDicrt2('e刁+尔本/2711rllNjIDOt=j-98根据结果可以得到，给出的迭代初值不一定会收敛于离它最近的实根，收敛速度也不一t定会慢

6、。初值不同所得到的收敛值也不同。例如，在本题中更改初值为x0=4时，所得到的根是3,迭代次数为4次。五、完成题目的体会与收获通过自己编程实现牛顿迭代法，不仅让我对牛顿迭代法有了更深刻的了解，同时也锻炼了我编程解决数学问题的能力。原本上课时不清晰的思路被理清了，观察计算结果之后，还对牛顿迭代法的规律和用法更加明了。希望以后能多有这样的实践作业。六、附录functionroot,t=NewtonRoot2(f,a)%f是非线性函数%a为初值%eps为根的精度%root为求出的函数零点%t为迭代次数eps=5.0e-6;t=0;f1=subs(sym(f),findsym(sym(f),a);fun

7、=diff(sym(f);fa=subs(sym(f),findsym(sym(f),a);dfa=subs(sym(fun),findsym(sym(fun),a);root=a-fa/dfa;tol=abs(root-a);while(tol>eps)t=t+1;r1=root;fx=subs(sym(f),findsym(sym(f),r1);dfx=subs(sym(fun),findsym(sym(fun),r1);root=r1-fx/dfx;tol=abs(root-r1);endend题目二：线性方程组求解有一平面机构如图所示，该机构共有13条梁（图中标号的线段）由8个较

8、接点（图中标号的圈）联结在一起。上述结构的1号钱接点完全固定，8号较接点竖立方向固定，并在2号、5号和6号较接点，分别有如图所示的10吨、15吨和20吨的负载，在静平衡的条件下，任何一个较接点上水平和竖立方向受力都是平衡的，以此计算每个梁的受力情况。令a=1/,假设f为各个梁上的受力，例如对2号钱接点有：f2=f6、f3=1。对3号钱接点有：otf,=af5+f4>0（1+otf5+f3=0对4号钱接点有：f4=f8、f7=0对5号钱接点有：Otfg+0tf5+f7=15、Ctf5+f6=afg+f10对6号皎接点有：10=f13、fn=20对7号钱接点有：af12=afg+f8、af1

9、2+ctfg+降=0对8号钱接点有：6后+%=0一、摘要对于实际的工程问题，很多问题归结为线性方程组的求解。本实验通过实际题目掌握求解线性方程组的数值解法，这里采用雅克比迭代法，如不收敛，再采用高斯列主元消去法。、数学原理1、雅克比迭代法设有一个n元线性方程组aiiXi-812X2一-am%a2lXi822X2一a2nXn=b2an1x1an2X2''annxnbn它的矩阵形式为AX=B,如果A=a)n约非奇异,且a.#0,i=1,2，,n。由上式可以得到Xi1，、(b-aijXj)(i=1,2,.,n)aiij4.jv而其相应的迭代公式为n1(k)、x=一(bi一2”)(i=

10、1,2,.,n)aiija及把上式迭代公式称为Jacobi(雅克比)迭代。由于迭代存在收敛性,an所以把分量形式的迭代公式改写成矩阵形式。记y0a2ia1201in2n-an2<ann八则A=DL-U.方程组Ax=b改写成x=D(LU)xDb与其相应的矩阵形式的迭代公式为X(k1)=D九LU)XkDb也可以简单地记为X(k+=BjXk+fJ式中，Bj=d,(l比)；fj=口节，上两式也称为Jacobi迭代。同时称Bj为Jacobi迭代矩阵。2、高斯列主元消去法在消元过程进行到第k步时，写出其相应的增广矩阵，可以发现，此时第k个方程与后面的n-k个方程的地位并没有区别，因此选择第k列的元素

11、a(kk)(i=k,kY.,n)中绝对值最大的元素作为主元，即令二息,瞟如果这时候a(kk)=0,那么矩阵就奇异不可逆，方程的解也不确定，只有停止计算；否则，当"则其增广矩阵中交换第k行和第r行，即akjk)1a：k)(j=k,k1,.,n)b，使a；?成为主元，然后再按高斯消去法进行消元运算。上述这种消去法称为高斯列主元消去法三、程序设计把方程组整理为矩阵形式:G010a000000001<0、01000-100000000001000000000010a00-1-a0000000000000a010a0000150000a100内-100000000001000000f=0

12、0001000-100000000000000a01a00000000000100-1000000000001002000000001a00-a00,100000000000a10°J本题我先采用了雅克比迭代法进行计算，所得结果发散，因此采用高斯列主元消去法计算。1、输入数据A和b,置det=1o2、对于k=1,2,.,n1作，按列选主元、交换两行、消元计算。»冀=b)x=7机3M320.MOO10.NOO-30.000014.H2120.QDOQ0-30.0DOO7.071125.000020.0000-35-366325.MOO3、置det=anndet。4、输出线性方

13、程组的解。四、结果分析和讨论得到结果，各个梁的受力情况分别为：-28.2843、20.000010.0000、-30.0000、14.1421、20.0000、0、-30.0000、7.071。25.0000、20.0000、-35.3553、25.0000（单位：吨）由结果分析，高斯列主消元法能准确的计算出该线性方程组的解。五、完成题目的体会与收获在解决本道题目的时候，我受到了重重困难。刚开始我并未考虑使用迭代法的收敛条件，便使用雅克比迭代法进行计算，但在经过多次尝试后，才发现该方法不收敛，改用高斯列主消元法来计算。这让我吸取了深深的教训。在今后的学习中，注意每种方法的使用限制条件，收敛条件

14、等，真正学而会用，才能彻底掌握知识。六、附录高斯列主消元法:functionx=Gauss(A,b)n,m=size(A);det=1;x=zeros(n,1);fork=1:n-1max1=0;fori=k:nifabs(A(i,k)>max1max1=abs(A(i,k);r=i;endendifr>kz=A(k,:);A(k,:)=A(r,:);A(r,:)=z;z=b(k);b(k)=b(r);b(r)=z;det=-det;endfori=k+1:nm=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);fork=n:-1:1forj=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=