哈尔滨工程大学数值分析大作业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 =时给出结果并分析现象,当6

510ε-=?，迭代停止。

二 数学原理

对于方程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）对于3

10x x --=，按照上述数学原理，编制的程序如下

program newton implicit none

real :: x(0:50),fx(0:50),f1x(0:50)!分别为自变量x ，函数f(x)和一阶导数f1(x)

integer :: k

write(*,*) "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,'(I3,1x,f11.6)') k,x(k)

if(abs(x(k)-x(k-1))<1e-6) exit !终止迭代条件 end do

stop end

（2）对于32

943892940x x x +-+=，按照上述数学原理，编制的程序如下

program newton implicit none

real :: x(0:50),fx(0:50),f1x(0:50)!分别为自变量x ，函数f(x)和一阶导数f1(x)

integer :: k

write(*,*) "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,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，但收敛速度明显不同。当初始值x0充分接近于方程的单根时，可保证迭代序列快速收敛。在本例中，初始值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 *的新的近似值。在本例中，当初始值x 0=2时，在这个点的切线方程与x 轴的交点恰为方程的一个根x=-98,故迭代了两次就得到了结果。

五 完成题目的体会与收获

对于牛顿法，当初始值x 0充分接近于方程的单根时，可保证迭代序列快速收敛。当方程有不止一个根时，所得结果与初始值的选取有关 。

题目二：线性方程组求解

一摘要

对于实际的工程问题，很多问题归结为线性方程组的求解。本实验通过实际题目掌握求解线性方程组的数值解法，直接法或间接法。

有一平面机构如图所示，该机构共有13条梁（图中标号的线段）由8个铰接点（图中标号的圈）联结在一起。上述结构的1号铰接点完全固定，8号铰接点竖立方向固定，并在2号、5号和6号铰接点，分别有如图所示的10吨、15吨和20吨的负载，在静平衡的条件下，任何一个铰接点上水平和竖立方向受力

都是平衡的，以此计算每个梁的受力情况。

10

15

20

令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 21212122221112

11

方法说明（以4阶为例）：

8

第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 none

integer,parameter :: n=13 ！输入量：系数阵的维数

dimension :: a(n,n),b(n),x(n)

double precision a,b,x

real,parameter :: m=0.7071 !m代表α= 1/

integer :: i,j

logical l

data((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(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 if

stop

end

subroutine 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 j=k,n

if (abs(a(i,j))>d) then

d=abs(a(i,j))

js(k)=j

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

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)

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 do

100 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