哈尔滨工程大学数值分析大作业附fortran程序.docx

1、哈尔滨工程大学数值分析大作业附fortran程序B班大作业要求：1使用统一封皮；2.上交大作业内容包含：一摘要二数学原理三程序设计（必须对输入变量、输出变量进行说明；编程无语言要求，但程序要求通过）四结果分析和讨论五完成题目的体会与收获3提交大作业的时间：本学期最后一次课，或考前答疑；过期不计入成绩；4提交方式：打印版一份；或手写大作业，但必须使用A4纸。5.撰写的程序需打印出来作为附录。纟爾滅N燈2 4镌程设计课程名称:目：学 号:姓 名:完成时间:题目一：非线性方程求根一摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。本实验通过使用常 用的求解方法二分法和Newton

2、法及改进的Newton法处理几个题目，分析并总结不同方法处理问 题的优缺点。观察迭代次数，收敛速度及初值选取对迭代的影响。用Newton法计算下列方程(1)疋7-1 = 0 ,初值分别为列 i ,入)=0.45, % =0.65；(2)x3+94x2-389x +294 = 0其三个根分别为1,3,-98。当选择初值.v0 = 2时给出结果并分析现 象，当 = 5xl(r迭代停止。二数学原理对于方程f(x)=0,如果f(x)是线性函数，则它的求根是很容易的。牛顿迭代法实质上是一种线 性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。设已知方程f(x)=O有近似根Xk(

3、假定F(xQHO),将函数f(x)在点Xk进行泰勒展开，有f(x) 2 f(xk)+f(xk)(x-xk)+ 于是方程f(x)=O可近似的表示为f(xk)+f(xk)(x-xk)=0这是个线性方程，记其根为Xk+1,贝心紳的计算公式为k=09lv这就是牛顿迭代法或简称牛顿法。三程序设计(本程序由Fortran语言编制)(1)对于/-尤-1 = ，按照上述数学原理，编制的程序如下program newtonimplicit nonereal: x(0:50),fx(0:50),flx(0:50)!分别为自变量 x,函数 f(x)和一阶导数 fl(x)integer : kwrite(*,*) ,

4、fx(0)=tfread(*,*) x(0)输入变量：初始值x(0)open(10,flle=,l.txt,)do k=l,50,lfx(k)=x(k-l)*3-x(k.l)-lflx(k)=3*x(k-l)*2-lx(k)=x(k-l)-fx(k)/flx(k)!牛顿法write(铝(I3,lxH6门 k,x(k) !输出变量：迭代次数k及x的值wHte(lV(I3Jx,fn.6)J k，x(k)if(abs(x(k)-x(k-l) do k=b50Jfx(k)=x(k-l)*3+94*x(k-l)*2-389*x(k-l)+294flx(k)=3*x(k-l)*2+188*x(k-l)-3

5、89x(k)=x(k-l)-fx(k)/flx(k)!牛顿法writeLUai)(十 口)kx(k) !输出变量：迭代次数k及x的值k/X(k)if(abs(x(k)-x(k-l)5e-6)exit !终止迭代条件end doend四结果分析和讨论(Q对于方程尤7-1 = 0,当初值分别为x0 = l,列=0.45,兀=0.65时，所得结果如下kXk初始值1初始值2初始值3()xo10.450.651XI1.500000-3.0121025.7915912X21.347826-2.0465173.9098533X31.325200-1.3958492.6869634X41.324718-0.9

6、162361.9264205X51.324718-0.3545281.5097046X6-1.4622501.3501837X7-0.9701851.3253028X8-0.4531211.3247189X9-2.1193701.32471810X10-1.44601211xn-0.95718212X12-0.43116813X13-1.89852714X14-1.29275915X15-0.82741716X16-0.12613717X17-1.04590918X18-0.56460119X19-14.65421020X20-9.78310721X21-6.54137022X22-4.387

7、30123X23-2.95878924X24-2.01102125X25-1.37128326X26-0.89570027X27 0.31076928X28-1.32340729X29-0.85459830X30-0.20847031X31-1.12909032X32-0.66518233X331.25644434X341.32950635X351.32473936X361.32471837X371.324718结果分析与讨论：从计算结果可以看到，当取的初始值不同时，虽然均得到了近似解 x*=1.324718,但收敛速度明显不同。当初始值xo充分接近于方程的单根时，可保证迭代序列快速 收敛。在

8、本例中，初始值1、0.65和045距方程的单根越来越远，故收敛速度越来越慢。(2)对于方程?+94?-389x+ 294 = 0,当初始值xo=2时计算结果如下kXk初始值0xo21XI-98.0000002X2-98.000000结果分析与讨论：牛顿法有明显的几何解释，方程f(x)=O的根X*可解释为曲线y=f(x)与X轴的 交点的横坐标。设xk是根的某个近似值，过曲线y=f(x)横坐标为xk的点Pk引曲线y=f(x)的切 线，并将该切线与x轴的交点坐标xk+i作为x*的新的近似值。在本例中，当初始值xo=2时，在这 个点的切线方程与x轴的交点恰为方程的一个根x=9&故迭代了两次就得到了结果

9、。五完成题目的体会与收获对于牛顿法，当初始值xo充分接近于方程的单根时，可保证迭代序列快速收敛。当方程有不 止一个根时，所得结果与初始值的选取有关。题目二：线性方程组求解一摘要对于实际的工程问题，很多问题归结为线性方程组的求解。本实验通过实际题目掌握求解线性 方程组的数值解法，直接法或间接法。有一平面机构如图所示，该机构共有13条梁（图中标号的线段）由8个较接点（图中标号的 圈）联结在一起。上述结构的1号较接点完全固定，8号较接点竖立方向固定，并在2号、5号和 6号较接点，分别有如图所示的吨、15吨和20吨的负载，在静平衡的条件下，任何一个较接点 上水平和竖立方向受力都是平衡的，以此计算每个梁

10、的受力情况。101520令a = U迈,假设/为各个梁上的受力，例如对8号枝接点有+d) thend=abs(a(ij)js(k)=jis=iend ifend doend do !把行绝对值最大的元素换到主元位置if (d+1.0=1.0) then1=0else !最大元素为0无解if(js(k)/=k) thendo i=l,nt=a(i9k)a(l,k)=a(IJs(k)a(ijs(k)=tend do !最大元素不在K行，K行end ifif(is/=k) thendo j=k,nt=a(kj)a(kj)=a(isj)a(isj)=t 咬换到K列end dot=b(k)b(k)=b(

11、is)b(is)=tend if !最大元素在主对角线上end if !消去if (l=0) thenwrlte(*,100)returnend ifdo j=k+l,na(kj)=a(kj)/a(k,k)end dob(k)=b(k)/a(k,k) !求三角矩阵do i=k+l,ndo j=k+l,na(ij)=a(ij)-a(i,k)*a(kj)end dob(i)=b(i)-a(i,k)*b(k)end doend doif (abs(a(n,n)+1.0=1.0) then1=0write(*400)returnend ifx(n)=b(n)/a(nji)do i=n-lj,-lt=0

12、.0do j=i+l,nt=t+a(ij)*x(j)end dox(i)=b(i)-tend do100 format(lx/fair)js(n)=ndo k=njlif (js(k)/=k) thent=x(k)x(k)=x(js(k)x(js(k)=tend ifend doreturnend四结果分析和讨论由程序计算的各个杆的受力如下:flf2f3f4fsf7-28.2820.0010.00-30.0014.1420.000.00f8f9fiofnfl2fl3-30.007.0725.0020.00 35.3625.00结果分析与讨论：列方程时假定各杆均受拉力，得到的结果有负值，说明力与

13、假设方向相反。 五完成题目的体会与收获高斯消去法虽然能够执行，但随便用akk（k）作为主元，有时会扩大误差，导致结果不可靠，甚 至严重失真，而高斯列主元消去法则不会有这种情况发生。最初采用的是高斯赛德尔迭代法解此矩阵，但是发现很难得到收敛解。因为必须满足迭代条件才可以，而找到满足条件的迭代矩阵非常困难，故最终选取了没有收敛限制的直接法解此矩阵。附录题目一程序：(1)program newtonimplicit nonereal: x(0:50)亦(O:5O),fBc(O:5O)!分别为自变量 x,函数 f(x)和一阶导数 fl(x)integer: kwrite(*?) ,lx(0)=Mrea

14、d(*/) x(0) !输入变量：初始值刈0)open(10,file=,l.txt,)do k=l/50/lfx(k)=x(k-l)3-x(k-l)-l!牛顿法flx(k)=3#x(k-l)*t2-lx( k)=x(k-l)-fx(k)/flx(k)write(叮(13幼1久6门 k,x(k)!输出变量：迭代次数k及x的值!终止迭代条件write(V(l3JxjM6) k,x(k)if(abs(x(k)-x(k-l)le-6) exit end do stopend(2)program newt onimplicit nonereal: x(0:50)x(0:50),flx(0:50)!分别

15、为自变量 x,函数 f(x)和一阶导数 fl(x)integer: kwrite(*/) ,x(0)=,1readC?) x(0) !输入变量：初始值刈0)open(10,file=,l.txt,)do k=l/50/lfx(k)=x(k-l)#*3+94*x(k-l)t*2-389#x(k-l)+294flx(k)=3#x(k-l)*t2+188*x(k-l)-389x(k)=x(kl)-fx(k)/flx(k) !牛顿法!输出变量：迭代次数k及）（的值!终止迭代条件write(MV(l3J)(Jll6) k,x(k)if(abs(x(k)-x(k-l)0.,0./1/0/0/-1/&0,0

16、,0000,00,0,00,&0,0,0/)j)m0jm0 &0,0,0oood0,0,0,0mi /b=0. !输入量：方程的右边项b(3)=10.b(9)=15.b(ll)=20.write(叮(13(13(f52x)r) (a(ij)沪 1山3)启 1J3) !输出矩阵call agausfaxJJs)if (l/=0) thenwri9(3(5f82,/)r)(:) !输出结果 end ifstop end subroutine agaus(a,bnx,js)dimension a(n,n),x(n),b(n),js(n)double precision a,b,x,t1=1 !逻辑变

17、fi：do k=l,n-ld=0.0do i=k,ndo j=k,nif (abs(a(ij)d) thend=abs(a(ijjs(k)=jis=iend ifend doend do !把行绝对值最大的元素换到主元位置if (d+1.0=1.0) then1=0else !最大元素为0无解if(js(k)/=k) thent=a(i,k)a(i/k)=a(ijs(k)a(i,js(k)=tend do !最大元素不在K行，K行end iflf(is/=k) thendot=a(k,j)a(k,j)=a(is,j)a(isj)=t 咬换到K列end dot=b(k)b(k)=b(is)b(i

18、s)=tend if !最大元素在主对角线上end if !消去if (l=0) thenwrit 屮 JOO)returnend ifdo j=k+l,na(kj)=a(kj)/a(k,k)end dob(k)=b(k)/a(k,k) !求三角矩阵do j=k+l/na(i/j)=a(M)-a(i/kra(kJ)end dob(i)=b(i)a(i* 广 b(k)end doend doif (abs(a(n,n)+1.0=1.0) then1=0write(*,100)returnend ifx(n)=b(n)/a(n/n)do 1=01,1,1t=0.0do j=i+l/nend dox(i)=b(i)-tend do100 format(lx/faila)js(n)=ndo k=n/l/-lif (js(k)/=k) thent=x(k)x(js(k)=tend if end do returnend精品文档word文档可以编借！谢谢下载!