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

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

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

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

B班大作业要求：

1•使用统一封皮；

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

一摘要

二数学原理

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

四结果分析和讨论

五完成题目的体会与收获

3•提交大作业的时间：

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

4•提交方式：

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

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

"纟爾滅N燈24

镌程设计

课程名称:

目：

学号:

姓名:

完成时间:

题目一：

非线性方程求根

一摘要

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

本实验通过使用常用的求解方法二分法和Newton法及改进的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（假定F（xQHO）,将函数f（x）在点Xk进行泰勒展开，有

f（x）2f（xk）+f（xk）（x-xk）+…

于是方程f（x）=O可近似的表示为

f（xk）+f（xk）（x-xk）=0

这是个线性方程，记其根为Xk+1,贝心紳的计算公式为

k=09l^v

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

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

（1）对于/-尤-1=°，按照上述数学原理，编制的程序如下

programnewton

implicitnone

real:

:

x（0:

50）,fx（0:

50）,flx（0:

50）!

分别为自变量x,函数f（x）和一阶导数fl（x）

integer:

:

k

write（*,*）,fx（0）=tf

read（*,*）x（0）［输入变量：

初始值x（0）

open（10,flle=,l.txt,）

dok=l,50,l

fx（k）=x（k-l）**3-x（k.l）-l

flx（k）=3*x（k-l）**2-l

x（k）=x（k-l）-fx（k）/flx（k）

!

牛顿法

write（铝（I3,lx』H・6门k,x（k）!

输出变量：

迭代次数k及x的值

wHte（l

if（abs（x（k）-x（k-l））

终止迭代条件

enddo

end

（2）对于?

+94x2-389a+294=0,按照上述数学原理，编制的程序如下

programnewtonimplicitnonereal:

:

x（0:

50bfx（0:

50）,flx（0:

50）!

分别为自变量x,函数f（x）和一阶导数fl（x）integer:

:

kwritef*/）'^（0）="

read（*/）x（0）!

输入变量：

初始值x（0）

open（10,file=ll.txt>）dok=b50J

fx（k）=x（k-l）**3+94*x（k-l）**2-389*x（k-l）+294

flx（k）=3*x（k-l）**2+188*x（k-l）-389

x（k）=x（k-l）-fx（k）/flx（k）

!

牛顿法

writeL'Uai）（十口⑹'）

kx（k）!

输出变量：

迭代次数k及x的值

k/X（k）

if（abs（x（k）-x（k-l））<5e-6）

exit!

终止迭代条件

enddo

end

四结果分析和讨论

（Q对于方程尤‘7-1=0,当初值分别为x0=l,列=0.45,兀=0.65时，所得结果如下

k

Xk

初始值1

初始值2

初始值3

（）

xo

1

0.45

0.65

1

XI

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

xn

-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,但收敛速度明显不同。

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

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

（2）对于方程?

+94?

-389x+294=0,当初始值xo=2时计算结果如下

k

Xk

初始值

0

xo

2

1

XI

-98.000000

2

X2

-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&故迭代了两次就得到了结果。

五完成题目的体会与收获

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

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

题目二：

线性方程组求解

一摘要

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

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

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

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

101520

令a=U迈,假设/为各个梁上的受力，例如对8号枝接点有

+<12=°

对5号较接点，则有

功5+人=劝+九）

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

二数学原理

高斯列主元消去法:

«12

…«ln

al\

■■•

…«2n

Ai

%

•…a州

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

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

**

0*

0*

n*

第2步消元——在增广矩阵（A,b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其

所在行与第二行交换，再对（A,b）做初等行变换使原方程组转化为^

*

*

0

0

*

*

*

*

*

*

*

54.

*

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

按"9X39X29X1的顺序回代求解出方程组的解。

针对各个较接点列方程:

对2号枝接点有

厶“°

对3号较接点有

砒一九一必=0

砒+人+必=°

对4号较接点有

对5号较接点有

扎_af®_仏=0

对6号较接点有

凡-九=o

对7号较接点有

对8号较接点有

九"0

£+好厂对2=0

%+人+砒2=0

«/12+/13=0

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

programmain

implicitnone

integer,parameter:

:

n=13!

输入量：

系数阵的维数

real:

:

js（n）

dimension:

:

a（n,n）9b（n）^（n）

doubleprecisiona,b,x

real,parameter:

:

m=0.7071!

m代表a=1/JT

integer:

:

ij

logical1

data（（a（ij）J=l,13）,i=l,13）/m,（M・,0・,m,0・,0・,（M）・,0・,0・,0・,0・,&!

输入量：

方程的系数阵

0・」・,0.,0・,0・9・1・,0・,0・,0・,（）・,0・,0・,0・，&

0・,0・,1・,0・,0・,0・,0・,0・90・,0・,0・,0・,0・,&

0・,0・,0•丄,0・,0“0・,・1・,0・,0・,0・,0・,0・，&

m,O・,O・,・l・,・m,O・,O・,O・,O・,O・,O・,O・,O・，&

0・f』・,0・,mJ』・,0.,・m,・l・,0・,0・,0•，&

0・,O・QQ,O・,0・J40・,0.f0•，&

0・g,0・,0・,0.,（M・,m,0・,0・,・in,0•，&0・,0・,0・,0・,m，0・,l・,0・,ni，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・90・,0・,1・,0・,0•，&

O・,O・,O・,O・,O・,O・,O・Q,m,O・丄,m,0•，&

O.,OmOmO.,O.,O.,O.,O.,O.,O.,0.411,1./

b=0.!

输入量：

方程的右边项

b⑶=10.

b（9）=15.

b（11）=20.

write（*,H（13（13（f5.2,lx）/）尸）（（a（iJ）J=l,13）4=1,13）!

输出矩阵

callagaus（a,b,n9xJJs）

if（l/=0）then

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

）!

输出结果

endif

stop

endsubroutineagaus（a9b,n,x,ljs）

dimensiona（n,n）以（n）,b（n）Js（n）

doubleprecisiona,b,x,t

1=1!

逻辑变量

dok=l,n-l

d=0.0

doi=k,n

doj=k,n

if（abs（a（ij））>d）then

d=abs（a（ij））

js（k）=j

is=i

endif

enddo

enddo!

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

if（d+1.0==1.0）then

1=0

else!

最大元素为0无解

if（js（k）/=k）then

doi=l,n

t=a（i9k）

a（l,k）=a（IJs（k））

a（ijs（k））=t

enddo!

最大元素不在K行，K行

endif

if（is/=k）then

doj=k,n

t=a（kj）

a（kj）=a（isj）

a（isj）=t咬换到K列

enddo

t=b（k）

b（k）=b（is）

b（is）=t

endif!

最大元素在主对角线上

endif!

消去

if（l==0）then

wrlte（*,100）

return

endif

doj=k+l,n

a（kj）=a（kj）/a（k,k）

enddo

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

求三角矩阵

doi=k+l,n

doj=k+l,n

a（ij）=a（ij）-a（i,k）*a（kj）

enddo

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

enddo

enddo

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

1=0

write（*400）

return

endif

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

doi=n-lj,-l

t=0.0

doj=i+l,n

t=t+a（ij）*x（j）

enddo

x（i）=b（i）-t

enddo

100format（lx/fair）

js（n）=n

dok=nj

if（js（k）/=k）then

t=x（k）

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

x（js（k））=t

endif

enddo

return

end

四结果分析和讨论

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

fl

f2

f3

f4

fs

f7

-28.28

20.00

10.00

-30.00

14.14

20.00

0.00

f8

f9

fio

fn

fl2

fl3

-30.00

7.07

25.00

20.00

・35.36

25.00

结果分析与讨论：

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

五完成题目的体会与收获

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

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

因为必须满足迭代条

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

附录

题目一程序：

（1）

programnewton

implicitnone

real:

:

x（0:

50）亦（O:

5O）,fBc（O:

5O）!

分别为自变量x,函数f（x）和一阶导数fl（x）

integer:

:

k

write（*?

）,lx（0）=M

read（*/）x（0）!

输入变量：

初始值刈0）

open（10,file=,l.txt,）

dok=l/50/l

fx（k）=x（k-l）^3-x（k-l）-l

!

牛顿法

flx（k）=3#x（k-l）*t2-l

x（k）=x（k-l）-fx（k）/flx（k）

write（叮（13』幼1久・6门k,x（k）

!

输出变量：

迭代次数k及x的值

!

终止迭代条件

write（«V（l3JxjM・6）'）k,x（k）

if（abs（x（k）-x（k-l））

end

（2）

programnewton

implicitnone

real:

:

x（0:

50）」x（0:

50）,flx（0:

50）!

分别为自变量x,函数f（x）和一阶导数fl（x）

integer:

:

k

write（*/）,,x（0）=,1

readC?

）x（0）!

输入变量：

初始值刈0）

open（10,file=,l.txt,）

dok=l/50/l

fx（k）=x（k-l）#*3+94*x（k-l）t*2-389#x（k-l）+294

flx（k）=3#x（k-l）*t2+188*x（k-l）-389

x（k）=x（kl）-fx（k）/flx（k）!

牛顿法

!

输出变量：

迭代次数k及）（的值

!

终止迭代条件

write（MV（l3J）（Jll・6）'）k,x（k）

if（abs（x（k）-x（k-l））<5e-6）exitenddostop

end

题目二程序：

programmain

implicitnone

integecparameter:

:

n=13!

输入：

ft：

系数阵的维数

real:

:

js（n）

dimension:

:

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

doubleprecisiona,b,xreal.parameter:

:

m=0.7071!

m代表a=

integer:

:

ij

logicalI

data（（a（i/j）/j=l/13）/i=l/13）/口0・丄0・皿040・曲・0曲・,&!

输入量：

方程的系数阵

0•丄』・0,0・,丄,0・,0・/）・,0・,0・,0・,0・・&

0・OJ・Od0・X）・X）・,0・,0・,0・OO,&

0・O/M・d0・X）・4X）・,0・,0・,0・,0・,&

0・・0・』・Oml・』・XVmO・X）・X）・,&

0./0./0./0.,0./0.,1/0./0./0./0./0./0./&

0・,0.j）・,0・,0・,0・Oj・m0・O,・rn/）・,&0・/）・』・/）・m（M・X）・m0・OX）・O,&

0.,0.,0./0./0./0./0>0.,0./1/0/0/-1/&

0・,0・,0・000』・,0・0,0・」・,0・0,&

0・,0・,0・』・/）・j）・』・』・m0・j・m0"&

0・,0・,0・oood0・,0・,0・,0・mi・/

b=0.!

输入量：

方程的右边项

b（3）=10.

b（9）=15.

b（ll）=20.

write（叮（13（13（f5・2』x"）r）（（a（ij）沪1山3）启1J3）!

输出矩阵

callagausfa^^xJJs）

if（l/=0）then

wri9（^（3（5f8・2,/）r））（（:

）!

输出结果endif

stopendsubroutineagaus（a,bnx」,js）

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

doubleprecisiona,b,x,t

1=1!

逻辑变fi：

dok=l,n-l

d=0.0

doi=k,n

doj=k,n

if（abs（a（ij））>d）then

d=abs（a（ij»

js（k）=j

is=i

endif

enddo

enddo!

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

if（d+1.0==1.0）then

1=0

else!

最大元素为0无解

if（js（k）/=k）then

t=a（i,k）

a（i/k）=a（ijs（k））

a（i,js（k））=t

enddo!

最大元素不在K行，K行

endif

lf（is/=k）then

do

t=a（k,j）

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

a（isj）=t咬换到K列

enddo

t=b（k）

b（k）=b（is）

b（is）=t

endif!

最大元素在主对角线上

endif!

消去

if（l==0）then

writ屮JOO）

return

endif

doj=k+l,n

a（kj）=a（kj）/a（k,k）

enddo

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

求三角矩阵

doj=k+l/n

a（i/j）=a（M）-a（i/kra（kJ）

enddo

b（i）=b（i）・a（i*广b（k）

enddo

enddo

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

1=0

write（*,100）

return

endif

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

do1=01,1,1

t=0.0

doj=i+l/n

enddo

x（i）=b（i）-t

enddo

100format（lx/faila）

js（n）=n

dok=n/l/-l

if（js（k）/=k）then

t=x（k）

x（js（k））=t

endifenddoreturn

end精品文档word文档可以编借！

谢谢下载!