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

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

《哈尔滨工程大学数值分析大作业附fortran程序.docx》由会员分享，可在线阅读，更多相关《哈尔滨工程大学数值分析大作业附fortran程序.docx（22页珍藏版）》请在冰点文库上搜索。

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

B班大作业要求：

1.使用统一封皮；

2.上交大作业内容包含：

一摘要

二数学原理

三程序设计（必须对输入变量、输出变量进行说明；编程无语言要求，但程序要求通过）

四结果分析和讨论

五完成题目的体会与收获

3.提交大作业的时间：

本学期最后一次课，或考前答疑；过期不计入成绩；

4.提交方式：

打印版一份；或手写大作业，但必须使用A4纸。

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

课程设计

课程名称：

设计题目：

学号：

姓名：

完成时间：

题目一：

非线性方程求根

一摘要

非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。

本实验通过使用常用的求解方法二分法和Newton法及改进的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的计算公式为

，k=0,1,2,…

这就是牛顿迭代法或简称牛顿法。

三程序设计（本程序由Fortran语言编制）

（1）对于

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

programnewton

implicitnone

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'）

dok=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!

终止迭代条件

enddo

stop

end

（2）对于

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

programnewton

implicitnone

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'）

dok=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!

终止迭代条件

enddo

stop

end

四结果分析和讨论

（1）对于方程

，当初值分别为

，

，

时，所得结果如下

k

xk

初始值1

初始值2

初始值3

0

x0

1

0.45

0.65

1

x1

1.500000

-3.012102

5.791591

2

x2

1.347826

-2.046517

3.909853

3

x3

1.325200

-1.395849

2.686963

4

x4

1.324718

-0.916236

1.926420

5

x5

1.324718

-0.354528

1.509704

6

x6

-1.462250

1.350183

7

x7

-0.970185

1.325302

8

x8

-0.453121

1.324718

9

x9

-2.119370

1.324718

10

x10

-1.446012

11

x11

-0.957182

12

x12

-0.431168

13

x13

-1.898527

14

x14

-1.292759

15

x15

-0.827417

16

x16

-0.126137

17

x17

-1.045909

18

x18

-0.564601

19

x19

-14.654210

20

x20

-9.783107

21

x21

-6.541370

22

x22

-4.387301

23

x23

-2.958789

24

x24

-2.011021

25

x25

-1.371283

26

x26

-0.895700

27

x27

-0.310769

28

x28

-1.323407

29

x29

-0.854598

30

x30

-0.208470

31

x31

-1.129090

32

x32

-0.665182

33

x33

1.256444

34

x34

1.329506

35

x35

1.324739

36

x36

1.324718

37

x37

1.324718

结果分析与讨论：

从计算结果可以看到，当取的初始值不同时，虽然均得到了近似解x*=1.324718，但收敛速度明显不同。

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

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

（2）对于方程

，当初始值x0=2时计算结果如下

k

xk

初始值

0

x0

2

1

x1

-98.000000

2

x2

-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充分接近于方程的单根时，可保证迭代序列快速收敛。

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

题目二：

线性方程组求解

一摘要

对于实际的工程问题，很多问题归结为线性方程组的求解。

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

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

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

48

135791112

261013

101520

令

，假设

为各个梁上的受力，例如对8号铰接点有

对5号铰接点，则有

针对各个铰接点，列出方程并求出各个梁上的受力。

二数学原理

高斯列主元消去法：

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

第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式：

第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：

第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：

按x4x3x2x1的顺序回代求解出方程组的解。

针对各个铰接点列方程：

对2号铰接点有

对3号铰接点有

对4号铰接点有

对5号铰接点有

对6号铰接点有

对7号铰接点有

对8号铰接点有

三程序设计（本程序由Fortran语言编制）

programmain

implicitnone

integer,parameter:

:

n=13！

输入量：

系数阵的维数

real:

:

js（n）

dimension:

:

a（n,n）,b（n）,x（n）

doubleprecisiona,b,x

real,parameter:

:

m=0.7071!

m代表α=

integer:

:

i,j

logicall

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（9）=15.

b（11）=20.

write（*,"（13（13（f5.2,1x）,/））"）（（a（i,j）,j=1,13）,i=1,13）！

输出矩阵

callagaus（a,b,n,x,l,js）

if（l/=0）then

write（*,"（3（5f8.2,/））"）x（:

）！

输出结果

endif

stop

end

subroutineagaus（a,b,n,x,l,js）

dimensiona（n,n）,x（n）,b（n）,js（n）

doubleprecisiona,b,x,t

l=1!

逻辑变量

dok=1,n-1

d=0.0

doi=k,n

doj=k,n

if（abs（a（i,j））>d）then

d=abs（a（i,j））

js（k）=j

is=i

endif

enddo

enddo!

把行绝对值最大的元素换到主元位置

if（d+1.0==1.0）then

l=0

else!

最大元素为0无解

if（js（k）/=k）then

doi=1,n

t=a（i,k）

a（i,k）=a（i,js（k））

a（i,js（k））=t

enddo!

最大元素不在K行，K行

endif

if（is/=k）then

doj=k,n

t=a（k,j）

a（k,j）=a（is,j）

a（is,j）=t!

交换到K列

enddo

t=b（k）

b（k）=b（is）

b（is）=t

endif!

最大元素在主对角线上

endif!

消去

if（l==0）then

write（*,100）

return

endif

doj=k+1,n

a（k,j）=a（k,j）/a（k,k）

enddo

b（k）=b（k）/a（k,k）!

求三角矩阵

doi=k+1,n

doj=k+1,n

a（i,j）=a（i,j）-a（i,k）*a（k,j）

enddo

b（i）=b（i）-a（i,k）*b（k）

enddo

enddo

if（abs（a（n,n））+1.0==1.0）then

l=0

write（*,100）

return

endif

x（n）=b（n）/a（n,n）

doi=n-1,1,-1

t=0.0

doj=i+1,n

t=t+a（i,j）*x（j）

enddo

x（i）=b（i）-t

enddo

100format（1x,'fail'）

js（n）=n

dok=n,1,-1

if（js（k）/=k）then

t=x（k）

x（k）=x（js（k））

x（js（k））=t

endif

enddo

return

end

四结果分析和讨论

由程序计算的各个杆的受力如下：

f1

f2

f3

f4

f5

f6

f7

-28.28

20.00

10.00

-30.00

14.14

20.00

0.00

f8

f9

f10

f11

f12

f13

-30.00

7.07

25.00

20.00

-35.36

25.00

结果分析与讨论：

列方程时假定各杆均受拉力，得到的结果有负值，说明力与假设方向相反。

五完成题目的体会与收获

高斯消去法虽然能够执行，但随便用akk（k）作为主元，有时会扩大误差，导致结果不可靠，甚至严重失真，而高斯列主元消去法则不会有这种情况发生。

最初采用的是高斯-赛德尔迭代法解此矩阵，但是发现很难得到收敛解。

因为必须满足迭代条件才可以，而找到满足条件的迭代矩阵非常困难，故最终选取了没有收敛限制的直接法解此矩阵。

附录

题目一程序：

（1）

programnewton

implicitnone

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'）

dok=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!

终止迭代条件

enddo

stop

end

（2）

programnewton

implicitnone

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'）

dok=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!

终止迭代条件

enddo

stop

end

题目二程序：

programmain

implicitnone

integer,parameter:

:

n=13！

输入量：

系数阵的维数

real:

:

js（n）

dimension:

:

a（n,n）,b（n）,x（n）

doubleprecisiona,b,x

real,parameter:

:

m=0.7071!

m代表α=

integer:

:

i,j

logicall

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（9）=15.

b（11）=20.

write（*,"（13（13（f5.2,1x）,/））"）（（a（i,j）,j=1,13）,i=1,13）！

输出矩阵

callagaus（a,b,n,x,l,js）

if（l/=0）then

write（*,"（3（5f8.2,/））"）x（:

）！

输出结果

endif

stop

end

subroutineagaus（a,b,n,x,l,js）

dimensiona（n,n）,x（n）,b（n）,js（n）

doubleprecisiona,b,x,t

l=1!

逻辑变量

dok=1,n-1

d=0.0

doi=k,n

doj=k,n

if（abs（a（i,j））>d）then

d=abs（a（i,j））

js（k）=j

is=i

endif

enddo

enddo!

把行绝对值最大的元素换到主元位置

if（d+1.0==1.0）then

l=0

else!

最大元素为0无解

if（js（k）/=k）then

doi=1,n

t=a（i,k）

a（i,k）=a（i,js（k））

a（i,js（k））=t

enddo!

最大元素不在K行，K行

endif

if（is/=k）then

doj=k,n

t=a（k,j）

a（k,j）=a（is,j）

a（is,j）=t!

交换到K列

enddo

t=b（k）

b（k）=b（is）

b（is）=t

endif!

最大元素在主对角线上

endif!

消去

if（l==0）then

write（*,100）

return

endif

doj=k+1,n

a（k,j）=a（k,j）/a（k,k）

enddo

b（k）=b（k）/a（k,k）!

求三角矩阵

doi=k+1,n

doj=k+1,n

a（i,j）=a（i,j）-a（i,k）*a（k,j）

enddo

b（i）=b（i）-a（i,k）*b（k）

enddo

enddo

if（abs（a（n,n））+1.0==1.0）then

l=0

write（*,100）

return

endif

x（n）=b（n）/a（n,n）

doi=n-1,1,-1

t=0.0

doj=i+1,n

t=t+a（i,j）*x（j）

enddo

x（i）=b（i）-t

enddo

100format（1x,'fail'）

js（n）=n

dok=n,1,-1

if（js（k）/=k）then

t=x（k）

x（k）=x（js（k））

x（js（k））=t

endif

enddo

return

end