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

B班大作业要求：1. 使用统一封皮；2. 上交大作业内容包含：一摘要二数学原理三程序设计（必须对输入变量、输出变量进行说明；编程无语言要求，但程序要求通过）四结果分析和讨论五完成题目的体会与收获3. 提交大作业的时间：本学期最后一次课，或考前答疑；过期不计入成绩；4. 提交方式：打印版一份；或手写大作业，但必须使用A4纸。

5. 撰写的程序需打印出来作为附录。

课程设计课程名称：设计题目：学号：姓名：完成时间：题目一：非线性方程求根 一 摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。

本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目，分析并总结不同方法处理问题的优缺点。

观察迭代次数，收敛速度及初值选取对迭代的影响。

用Newton 法计算下列方程（1） 310x x --= ， 初值分别为01x =，00.45x =，00.65x =；（2） 32943892940x x x +-+= 其三个根分别为1,3,98-。

当选择初值02x =时给出结果并分析现象,当6510ε-=⨯，迭代停止。

二 数学原理对于方程f(x)=0，如果f(x)是线性函数，则它的求根是很容易的。

牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。

设已知方程f(x)=0有近似根x k (假定k f'(x )0≠) ，将函数f(x)在点x k 进行泰勒展开，有k k k f(x)f(x )+f'(x )(x-x )+≈⋅⋅⋅于是方程f(x)=0可近似的表示为k k k f(x )+f'(x )(x-x )=0这是个线性方程，记其根为x k+1，则x k+1的计算公式为k+1k ()x =x -'()k k f x f x ，k=0,1,2,… 这就是牛顿迭代法或简称牛顿法。

三 程序设计（本程序由Fortran 语言编制）（1）对于310x x --=，按照上述数学原理，编制的程序如下program newton implicit 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,1fx(k)=x(k-1)**3-x(k-1)-1 f1x(k)=3*x(k-1)**2-1x(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 do stop end（2）对于32943892940x x x +-+=，按照上述数学原理，编制的程序如下program newton implicit 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,1fx(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)-389x(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 do stop end四 结果分析和讨论（1）对于方程310x x --=，当初值分别为01x =，00.45x =，00.65x =时，所得结果如下结果分析与讨论：从计算结果可以看到，当取的初始值不同时，虽然均得到了近似解x *=1.324718，但收敛速度明显不同。

当初始值x 0充分接近于方程的单根时，可保证迭代序列快速收敛。

在本例中，初始值1、0.65和0.45距方程的单根越来越远，故收敛速度越来越慢。

（2）对于方程32943892940x x x +-+=，当初始值x 0=2时计算结果如下结果分析与讨论：牛顿法有明显的几何解释，方程f(x)=0的根x*可解释为曲线y=f(x)与x轴的交点的横坐标。

设x k是根x*的某个近似值，过曲线y=f(x)上横坐标为x k的点P k引曲线y=f(x)的切线，并将该切线与x轴的交点坐标x k+1作为x*的新的近似值。

在本例中，当初始值x0=2时，在这个点的切线方程与x轴的交点恰为方程的一个根x=-98,故迭代了两次就得到了结果。

五完成题目的体会与收获对于牛顿法，当初始值x0充分接近于方程的单根时，可保证迭代序列快速收敛。

当方程有不止一个根时，所得结果与初始值的选取有关。

题目二：线性方程组求解一摘要对于实际的工程问题，很多问题归结为线性方程组的求解。

本实验通过实际题目掌握求解线性方程组的数值解法，直接法或间接法。

有一平面机构如图所示，该机构共有13条梁（图中标号的线段）由8个铰接点（图中标号的圈）联结在一起。

上述结构的1号铰接点完全固定，8号铰接点竖立方向固定，并在2号、5号和6号铰接点，分别有如图所示的10吨、15吨和20吨的负载，在静平衡的条件下，任何一个铰接点上水平和竖立方向受力都是平衡的，以此计算每个梁的受力情况。

101520令2/1=α，假设f 为各个梁上的受力，例如对8号铰接点有01213=+f f α对5号铰接点，则有10965f f f f +=+αα15975=++f f f αα针对各个铰接点，列出方程并求出各个梁上的受力。

二 数学原理 高斯列主元消去法：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n n n b b b x x x a a a a a a a a a 2121212222111211方法说明（以4阶为例）：第1步消元——在增广矩阵（A ，b ）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A ，b ）做初等行变换使原方程组转化为如下形式：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡*******0***0***0****4321x x x x第2步消元——在增广矩阵（A ，b ）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A ，b ）做初等行变换使原方程组转化为：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡******00**00***0****4321x x x x第3步消元——在增广矩阵（A ，b ）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A ，b ）做初等行变换使原方程组转化为：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡*****000**00***0****4321x x x x按x 4→ x 3→ x 2→ x 1 的顺序回代求解出方程组的解。

针对各个铰接点列方程： 对2号铰接点有260f f -= 310f =对3号铰接点有1450f f f αα--=1350f f f αα++=对4号铰接点有480f f -=70f =对5号铰接点有569100f f f f αα+--=579150f f f αα++-=对6号铰接点有10130f f -= 1120f =对7号铰接点有89120f f f αα+-=911120f f f αα++=对8号铰接点有12130f f α+=三 程序设计（本程序由Fortran 语言编制）program main implicit 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 代表α=1/integer :: i,j logical 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.,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.,-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) thenwrite(*,"(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,tl=1 !逻辑变量do k=1,n-1d=0.0do i=k,ndo j=k,nif (abs(a(i,j))>d) thend=abs(a(i,j))js(k)=jis=iend ifend doend do !把行绝对值最大的元素换到主元位置if (d+1.0==1.0) thenl=0else !最大元素为0无解if(js(k)/=k) thendo i=1,nt=a(i,k)a(i,k)=a(i,js(k))a(i,js(k))=tend do !最大元素不在K行，K行end ifif(is/=k) thendo j=k,nt=a(k,j)a(k,j)=a(is,j)a(is,j)=t !交换到K列end dot=b(k)b(k)=b(is)b(is)=tend if !最大元素在主对角线上end if !消去if (l==0) thenwrite(*,100)returnend ifdo j=k+1,na(k,j)=a(k,j)/a(k,k)end dob(k)=b(k)/a(k,k) !求三角矩阵do i=k+1,ndo j=k+1,na(i,j)=a(i,j)-a(i,k)*a(k,j)end dob(i)=b(i)-a(i,k)*b(k)end doend doif (abs(a(n,n))+1.0==1.0) thenl=0write(*,100)returnend ifx(n)=b(n)/a(n,n)do i=n-1,1,-1t=0.0do j=i+1,nt=t+a(i,j)*x(j)end dox(i)=b(i)-tend do100 format(1x,'fail')js(n)=ndo k=n,1,-1if (js(k)/=k) thent=x(k)x(k)=x(js(k))x(js(k))=tend ifend doreturnend四结果分析和讨论由程序计算的各个杆的受力如下：结果分析与讨论：列方程时假定各杆均受拉力，得到的结果有负值，说明力与假设方向相反。