数值分析报告大作业2014

1、word用演？免火孝课程设计课程名称：高等数值计算设计题目：数值计算B课程设计学号：某某：完成时间：2014年10月20日1 / 12word题目一：非线性方程求根用Newton法计算如下方程1x3X 10,初值分别为X01,X00.45,x。0.65.2x3 94x2 389x 294 0其三个根分别为13 98。当选择初值xo 2时给 出结果并分析现象，当 5 106,迭代停止。一、摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。本实验通过使用常用的求解方法二分法和 Newton法与改良的Newton法处理几 个题目，分析并总结不同方法处理问题的优缺点。 观察迭代次数

2、，收敛速度与初 值选取对迭代的影响。二、数学原理构造迭代函数的一条很重要的途径是，用近似方程来代替原方程去求根。因 此，如果能将非线性方程用线性方程来代替的话，求近似根问题就很容易解决， 而且十分方便。Newton法就是把非线性方程线性化的一种方法。在求解非线性方程f(x) 0时，它的困难在于f(x)是非线性函数，为克制这 一困难，考虑它的线性展开。设当前点为 xk,在xk处的Taylor展开式为f (x) f (xk) f (xk)(x xk)令f(x) 0,可以得到上式的近似方程f (xk) f (xk)(x xk) 02 / 12word设f(%)。，解其方程得到f (凡)Xki Xk(

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

4、计算所得Xk 1-Xk是否满足精度要求4、如不满足继续迭代运算，如满足如此输出所求结果四、结果分析和讨论1、第一题计算结果：首先得到函数在区间，2.5的图像，即可知函数f(x)与x轴有交点，也就是 说有根，并且从图中能够大致估算到根的位置。3 / 12wordi3m中用工*式"T才Arilxotrt =t =4»J 二 Be =聋rod I.有t - 2» 【not/】 HiHtaidCQan" " "X"Il、 rcio4 .1. Rd，(1)、取初值Xo 1时得到根值，迭代次数t=4次(2)、取初值Xo 0.45时得到根

5、值，迭代次数t=42次(3)、取初值Xo 0.65时得到根值，迭代次数t=8次根据结果可以分析得到，当使用牛顿迭 代法时，所选初始值对迭代速度迭代次数 有较大影响。当初始值Xo充分接近方程的单根 时，可保证迭代序列快速收敛，当初值选择不 当时会造成迭代次数大幅增加或不一定收敛。2、第二题计算结果：初值Xo 2时，得到根值r=-98,迭代次数” IrEtt - ffevt rnRiMrt 2( l葭 3+与 *-.？ 2为1次。皿1 0-S3根据结果可以得到，给出的迭代初值不一定会收敛于离它最近的实根，收敛速度也不一一定会慢。初值不同所得到的收敛值也不同。例如，在此题中更改初值为Xo 4时，所得

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

7、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;5 / 12wordtol=abs(root-r1);endend题目二：线性方程组求解有一平面机构如下列图，该机构共有 13条梁图中标号的线段由8个较接点图中标号的圈联结在一

8、起。上述结构的 1号较接点完全固定，8号较接点竖 立方向固定，并在2号、5号和6号较接点，分别有如下列图的10吨、15吨和 20吨的负载，在静平衡的条件下，任何一个较接点上水平和竖立方向受力都是 平衡的，以此计算每个梁的受力情况。6 / 12word101520令1/J2,假设f为各个梁上的受力，例如对2号钱接点有：f2 f6> f3 10对3号钱接点有：f1f5 f4、f1f5 f3 0对4号钱接点有：f4 f8、f7 0对 5 号钱接点有：f9f5 f7 15、f5 f6f9 f10对6号钱接点有：f10 f13、fn 20对 7 号钱接点有：f12f9 f8、f12f9 f11 0

9、对8号钱接点有：f12 f13 0一、摘要对于实际的工程问题，很多问题归结为线性方程组的求解。本实验通过实际 题目掌握求解线性方程组的数值解法， 这里采用雅克比迭代法，如不收敛，再采 用高斯列主元消去法。二、数学原理1、雅克比迭代法设有一个n元线性方程组a11X1 812X2 anXn 6a21X1822X2 a2nXn b20,i 1,2，,n。由上式an1X1 an2X2 annXn它的矩阵形式为AX B,如果A (aij)nn非奇异，且叫7 / 12可以得到xi 一(baiiwordaijxj) (i 1,2,.,n)11而其相应的迭代公式为x (bia4Xj(k) (i 1,2,.,

10、n)j 1j 1把上式迭代公式称为Jacobi(雅克比)迭代。由于迭代存在收敛性，所以把分量形式的迭代公式改写成矩阵形式。记即a2200a12.a1ra210,0.a2.an1.an2.0.0nnann如此A D L U .方程组Ax b改写成x D 1(LU)xD 1b与其相应的矩阵形式的迭代公式为x(k 1) D 1(L U )xk D 1b也可以简单地记为(k 1)kxBjX fj式中，Bj D 1(L U) ; fjd 1b，上两式也称为Jacobi迭代。同时称Bj为Jacobi迭代矩阵。2、高斯列主元消去法在消元过程进展到第k步时，写出其相应的增广矩阵，可以发现，此时第k个方程与后面

11、的n k个方程的地位并没有区别，因此选择第k列的元素 ai(kk)(ik,k 1,.,n)中绝对值最大的元素作为主元，即令a：max 鼠"k i n8 / 12word如果这时候a(kk)=0 ,那么矩阵就奇异不可逆，方程的解也不确定，只有停止 计算；否如此，当r k,如此其增广矩阵换第k行和第r行，即a；k)a：k)(j k,k 1,.,n)使a会成为主元，然后再按高斯消去法进展消元运算。上述这种消去法称为 高斯列主元消去法。三、程序设计把方程组整理为矩阵形式:000010 000000000001010001000010000000000000010000000000000000

12、0000000010000000 f00000101001010000000100100150 000020 000 1 00 10 000 0 1000 010000000000000000100000000000000000000000000此题我先采用了雅克比迭代法进展计算,所得结果发散，因此采用高斯列主元消去法计算1、输入数据A和b,置det=1o2、对于k 1,2,., n 1作，按列选主元、交换两行、消元计算3、置 det ann det。4、输出线性方程组的解。四、结果分析和讨论得到结果，各个梁的受力情况分别为：-28.2843、20.000010.0000、-30.0000、

13、14.1421、20.0000、0、-30.0000、7.0711、-51. 94320.QOOQ10. fflHOi-30-000Q 露 U2I 也 aooo7.07 I2n. MQQ7 口知25.QOOQ9 / 12word25.0000、20.0000、-35.3553、25.0000单位：吨由结果分析，高斯列主消元法能准确的计算出该线性方程组的解。五、完成题目的体会与收获在解决本道题目的时候，我受到了重重困难。刚开始我并未考虑使用迭代法 的收敛条件，便使用雅克比迭代法进展计算，但在经过屡次尝试后，才发现该方 法不收敛，改用高斯列主消元法来计算。这让我吸取了深深的教训。在今后的学 习中，

14、注意每种方法的使用限制条件，收敛条件等，真正学而会用，才能彻底掌 握知识。六、附录高斯列主消元法：function x= Gauss(A,b)n,m=size(A);det=1;x=zeros(n,1);for k=1:n-1max1=0;for i=k:nif abs(A(i,k)>max1max1=abs(A(i,k);r=i;endendif r>kz=A(k,:);A(k,:)=A(r,:);A(r,:)=z;z=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n10 / 12wordm=A(i,k)/A(k,k);for j=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);for k=n:-1:1for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);endx(k)=b