哈尔滨工建大学数值大作业2014-附fortran程序

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

2、进的Newton法处理几个题目，分析并总结不同方法处理问题的优缺点。观察迭代次数，收敛速度及初值选取对迭代的影响。用Newton法计算下列方程 （1） ， 初值分别为，； （2） 其三个根分别为。当选择初值时给出结果并分析现象,当，迭代停止。二 数学原理对于方程f(x)=0，如果f(x)是线性函数，则它的求根是很容易的。牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。设已知方程f(x)=0有近似根xk(假定) ，将函数f(x)在点xk进行泰勒展开，有于是方程f(x)=0可近似的表示为这是个线性方程，记其根为xk+1，则xk+1的计算公式为 ，

3、k=0,1,2,这就是牛顿迭代法或简称牛顿法。三 程序设计（本程序由Fortran语言编制）（1）对于，按照上述数学原理，编制的程序如下 program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数f1(x)integer : kwrite(*,*) "x(0)="read(*,*) x(0) !输入变量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3-x(k-1)-1 f1x(k)=3*x(k-

4、1)*2-1 x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法 write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k及x的值 write(10,'(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1)<1e-6) exit !终止迭代条件end dostopend（2）对于，按照上述数学原理，编制的程序如下program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数f1(x)intege

5、r : kwrite(*,*) "x(0)="read(*,*) x(0) !输入变量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3+94*x(k-1)*2-389*x(k-1)+294 f1x(k)=3*x(k-1)*2+188*x(k-1)-389 x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法 write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k及x的值 write(10,'(I3,1x,f11.6)') k,

6、x(k) if(abs(x(k)-x(k-1)<5e-6) exit !终止迭代条件end dostopend四 结果分析和讨论（1）对于方程，当初值分别为，时，所得结果如下kxk初始值1初始值2初始值30x010.450.651x11.500000-3.0121025.7915912x21.347826-2.0465173.9098533x31.325200-1.3958492.6869634x41.324718-0.9162361.9264205x51.324718-0.3545281.5097046x6-1.4622501.3501837x7-0.9701851.3253028x8

7、-0.4531211.3247189x9-2.1193701.32471810x10-1.44601211x11-0.95718212x12-0.43116813x13-1.89852714x14-1.29275915x15-0.82741716x16-0.12613717x17-1.04590918x18-0.56460119x19-14.65421020x20-9.78310721x21-6.54137022x22-4.38730123x23-2.95878924x24-2.01102125x25-1.37128326x26-0.89570027x27-0.31076928x28-1.32

8、340729x29-0.85459830x30-0.20847031x31-1.12909032x32-0.66518233x331.25644434x341.32950635x351.32473936x361.32471837x371.324718结果分析与讨论：从计算结果可以看到，当取的初始值不同时，虽然均得到了近似解x*=1.324718，但收敛速度明显不同。当初始值x0充分接近于方程的单根时，可保证迭代序列快速收敛。在本例中，初始值1、0.65和0.45距方程的单根越来越远，故收敛速度越来越慢。（2）对于方程，当初始值x0=2时计算结果如下kxk初始值0x021x1-98.000000

9、2x2-98.000000结果分析与讨论：牛顿法有明显的几何解释，方程f(x)=0的根x*可解释为曲线y=f(x)与x轴的交点的横坐标。设xk是根x*的某个近似值，过曲线y=f(x)上横坐标为xk的点Pk引曲线y=f(x)的切线，并将该切线与x轴的交点坐标xk+1作为x*的新的近似值。在本例中，当初始值x0=2时，在这个点的切线方程与x轴的交点恰为方程的一个根x=-98,故迭代了两次就得到了结果。五 完成题目的体会与收获 对于牛顿法，当初始值x0充分接近于方程的单根时，可保证迭代序列快速收敛。当方程有不止一个根时，所得结果与初始值的选取有关 。 题目二：线性方程组求解一 摘要对于实际的工程问题

10、，很多问题归结为线性方程组的求解。本实验通过实际题目掌握求解线性方程组的数值解法，直接法或间接法。有一平面机构如图所示，该机构共有13条梁（图中标号的线段）由8个铰接点（图中标号的圈）联结在一起。上述结构的1号铰接点完全固定，8号铰接点竖立方向固定，并在2号、5号和6号铰接点，分别有如图所示的10吨、15吨和20吨的负载，在静平衡的条件下，任何一个铰接点上水平和竖立方向受力都是平衡的，以此计算每个梁的受力情况。7865434813579111221261013101520 令，假设为各个梁上的受力，例如对8号铰接点有对5号铰接点，则有针对各个铰接点，列出方程并求出各个梁上的受力。二 数学原理高

11、斯列主元消去法：方法说明（以4阶为例）：第1步消元在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式：第2步消元在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：第3步消元在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：按x4 à x3à x2à x1 的顺序回代求解出方程组的解。针对各个铰接点列方程：对2号铰接点有

12、对3号铰接点有对4号铰接点有对5号铰接点有对6号铰接点有对7号铰接点有对8号铰接点有三 程序设计（本程序由Fortran语言编制）program mainimplicit noneinteger,parameter : n=13 ！输入量：系数阵的维数real : js(n)dimension : a(n,n),b(n),x(n)double precision a,b,xreal,parameter : m=0.7071 !m代表= integer : i,jlogical ldata(a(i,j),j=1,13),i=1,13) / m,0.,1.,0.,m,0.,0.,0.,0.,0.,

13、0.,0.,0., & ！输入量：方程的系数阵 0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0., & m,0.,0.,-1.,-m,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,0.,m,1.,0.,0.,-m,-1.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,&a

14、mp;0.,0.,0.,0.,0.,0.,0.,1.,m,0.,0.,-m,0.,& 0.,0.,0.,0.,m,0.,1.,0.,m,0.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,-1.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,m,0.,1.,m,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,m,1. / b=0. ！输入量：方程的右边项 b(3)=10. b(9)=15. b

15、(11)=20.write(*,"(13(13(f5.2,1x),/)") (a(i,j),j=1,13),i=1,13) ！输出矩阵call agaus(a,b,n,x,l,js) if (l/=0) then write(*,"(3(5f8.2,/)") x(:) ！输出结果end ifstopendsubroutine agaus(a,b,n,x,l,js) dimension a(n,n),x(n),b(n),js(n) double precision a,b,x,t l=1 !逻辑变量 do k=1,n-1 d=0.0 do i=k,n do

16、 j=k,n if (abs(a(i,j)>d) then d=abs(a(i,j) js(k)=j is=i end if end do end do !把行绝对值最大的元素换到主元位置 if (d+1.0=1.0) then l=0 else !最大元素为0无解 if(js(k)/=k) then do i=1,n t=a(i,k) a(i,k)=a(i,js(k) a(i,js(k)=t end do !最大元素不在K行，K行 end if if(is/=k) then do j=k,n t=a(k,j) a(k,j)=a(is,j) a(is,j)=t !交换到K列 end do

17、 t=b(k) b(k)=b(is) b(is)=t end if !最大元素在主对角线上 end if !消去 if (l=0) then write(*,100) return end if do j=k+1,n a(k,j)=a(k,j)/a(k,k) end do b(k)=b(k)/a(k,k) !求三角矩阵 do i=k+1,n do j=k+1,n a(i,j)=a(i,j)-a(i,k)*a(k,j) end do b(i)=b(i)-a(i,k)*b(k) end do end do if (abs(a(n,n)+1.0=1.0) then l=0 write(*,100)

18、return end if x(n)=b(n)/a(n,n) do i=n-1,1,-1 t=0.0 do j=i+1,n t=t+a(i,j)*x(j) end do x(i)=b(i)-t end do100 format(1x,'fail') js(n)=n do k=n,1,-1 if (js(k)/=k) then t=x(k) x(k)=x(js(k) x(js(k)=t end if end do return end四 结果分析和讨论由程序计算的各个杆的受力如下：f1f2f3f4f5f6f7-28.2820.0010.00-30.0014.1420.000.00

19、f8f9f10f11f12f13-30.007.0725.0020.00-35.3625.00结果分析与讨论：列方程时假定各杆均受拉力，得到的结果有负值，说明力与假设方向相反。五 完成题目的体会与收获高斯消去法虽然能够执行，但随便用akk(k)作为主元，有时会扩大误差，导致结果不可靠，甚至严重失真，而高斯列主元消去法则不会有这种情况发生。最初采用的是高斯-赛德尔迭代法解此矩阵，但是发现很难得到收敛解。因为必须满足迭代条件才可以，而找到满足条件的迭代矩阵非常困难，故最终选取了没有收敛限制的直接法解此矩阵。附录题目一程序：（1）program newtonimplicit nonereal : x

20、(0:50),fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数f1(x)integer : kwrite(*,*) "x(0)="read(*,*) x(0) !输入变量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3-x(k-1)-1 f1x(k)=3*x(k-1)*2-1 x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法 write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k及x的值 write(10,&#

21、39;(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1)<1e-6) exit !终止迭代条件end dostopend（2）program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数f1(x)integer : kwrite(*,*) "x(0)="read(*,*) x(0) !输入变量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3+94*x(

22、k-1)*2-389*x(k-1)+294 f1x(k)=3*x(k-1)*2+188*x(k-1)-389 x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法 write(*,'(I3,1x,f11.6)') k,x(k) !输出变量：迭代次数k及x的值 write(10,'(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1)<5e-6) exit !终止迭代条件end dostopend题目二程序：program mainimplicit noneinteger,parameter : n=13 ！输入量：系数阵的

23、维数real : js(n)dimension : a(n,n),b(n),x(n)double precision a,b,xreal,parameter : m=0.7071 !m代表= integer : i,jlogical ldata(a(i,j),j=1,13),i=1,13) / m,0.,1.,0.,m,0.,0.,0.,0.,0.,0.,0.,0., & ！输入量：方程的系数阵 0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., & 0.,

24、0.,0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0., & m,0.,0.,-1.,-m,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,0.,m,1.,0.,0.,-m,-1.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,&0.,0.,0.,0.,0.,0.,0.,1.,m,0.,0.,-m,0.,& 0.,0.,0.,0.,m,0.,1.,0.,m,0.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0

25、.,-1.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,m,0.,1.,m,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,m,1. / b=0. ！输入量：方程的右边项 b(3)=10. b(9)=15. b(11)=20.write(*,"(13(13(f5.2,1x),/)") (a(i,j),j=1,13),i=1,13) ！输出矩阵call agaus(a,b,n,x,l,js) if (l/=0) then write(*,"(3(5f8.2,/)") x(:) ！输出结果end ifstopendsubroutine agaus(a,b,n,x,l,js) dimension a(n,n),x(n),b(n),js(n) double precision a,b,x,t l=1 !逻辑变量 do k