哈尔滨工程大学数值分析大作业附fortran程序Word文档格式.docx
《哈尔滨工程大学数值分析大作业附fortran程序Word文档格式.docx》由会员分享，可在线阅读，更多相关《哈尔滨工程大学数值分析大作业附fortran程序Word文档格式.docx（18页珍藏版）》请在冰豆网上搜索。
:
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<
V（I3Jx,fn.6）Jk，x（k）
if（abs（x（k）-x（k-l））<
le-6）exit!
终止迭代条件
enddo
end
（2）对于?
+94x2-389a+294=0,按照上述数学原理，编制的程序如下
programnewtonimplicitnonereal:
50bfx（0:
分别为自变量x,函数f（x）和一阶导数fl（x）integer:
kwritef*/）'
^（0）="
read（*/）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
writeL'
Uai）（十口⑹'
）
kx（k）!
k/X（k）
5e-6）
exit!
（Q对于方程尤‘7-1=0,当初值分别为x0=l,列=0.45,兀=0.65时，所得结果如下
k
Xk
初始值1
初始值2
初始值3
（）
xo
1
0.45
0.65
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
-0.354528
1.509704
6
X6
-1.462250
1.350183
7
X7
-0.970185
1.325302
8
X8
-0.453121
9
X9
-2.119370
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
37
X37
结果分析与讨论：
从计算结果可以看到，当取的初始值不同时，虽然均得到了近似解x*=1.324718,但收敛速度明显不同。
当初始值xo充分接近于方程的单根时，可保证迭代序列快速收敛。
在本例中，初始值1、0.65和0・45距方程的单根越来越远，故收敛速度越来越慢。
（2）对于方程?
+94?
-389x+294=0,当初始值xo=2时计算结果如下
初始值
-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+人=劝+九）
针对各个较接点，列出方程并求出各个梁上的受力。
高斯列主元消去法:
«
…«
ln
al\
■■•
2n
Ai
%
•…a州
方法说明（以4阶为例）：
第1步消元——在增广矩阵（A,b）第一列中找到绝对值最大的元素，将其所在行与第一行交换,再对（A,b）做初等行变换使原方程组转化为如下形式：
**
0*
n*
第2步消元——在增广矩阵（A,b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其
所在行与第二行交换，再对（A,b）做初等行变换使原方程组转化为^
*
54.
第3步消元一在增广矩阵（A,b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A,b）做初等行变换使原方程组转化为：
按"
9X39X29X1的顺序回代求解出方程组的解。
针对各个较接点列方程:
对2号枝接点有
厶“°
对3号较接点有
砒一九一必=0
砒+人+必=°
对4号较接点有
对5号较接点有
扎_af®
_仏=0
对6号较接点有
凡-九=o
对7号较接点有
对8号较接点有
九"
£
+好厂对2=0
%+人+砒2=0
/12+/13=0
三程序设计（本程序由Fortran语言编制）
programmain
integer,parameter:
n=13!
输入量：
系数阵的维数
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
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
最大元素不在K行，K行
if（is/=k）then
t=a（kj）
a（kj）=a（isj）
a（isj）=t咬换到K列
t=b（k）
b（k）=b（is）
b（is）=t
endif!
最大元素在主对角线上
消去
if（l==0）then
wrlte（*,100）
return
doj=k+l,n
a（kj）=a（kj）/a（k,k）
b（k）=b（k）/a（k,k）!
求三角矩阵
doi=k+l,n
a（ij）=a（ij）-a（i,k）*a（kj）
b（i）=b（i）-a（i,k）*b（k）
if（abs（a（n,n））+1.0==1.0）then
write（*400）
x（n）=b（n）/a（nji）
doi=n-lj,-l
t=0.0
doj=i+l,n
t=t+a（ij）*x（j）
x（i）=b（i）-t
100format（lx/fair）
js（n）=n
dok=nj<
l
if（js（k）/=k）then
t=x（k）
x（k）=x（js（k））
x（js（k））=t
由程序计算的各个杆的受力如下:
fl
f2
f3
f4
fs
f7
-28.28
20.00
10.00
-30.00
14.14
0.00
f8
f9
fio
fn
fl2
fl3
7.07
25.00
・35.36
列方程时假定各杆均受拉力，得到的结果有负值，说明力与假设方向相反。
五完成题目的体会与收获
高斯消去法虽然能够执行，但随便用akk（k）作为主元，有时会扩大误差，导致结果不可靠，甚至严重失真，而高斯列主元消去法则不会有这种情况发生。
最初采用的是高斯•赛德尔迭代法解此矩阵，但是发现很难得到收敛解。
因为必须满足迭代条
件才可以，而找到满足条件的迭代矩阵非常困难，故最终选取了没有收敛限制的直接法解此矩阵。
附录
题目一程序：
（1）
50）亦（O:
5O）,fBc（O:
5O）!
write（*?
）,lx（0）=M
初始值刈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）
write（«
V（l3JxjM・6）'
）k,x（k）
le-6）exitenddostop
（2）
programnewton
50）」x（0:
write（*/）,,x（0）=,1
readC?
）x（0）!
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）'
5e-6）exitenddostop
题目二程序：
integecparameter:
输入：
ft：
a（n,n）,b（n）,x（n）
doubleprecisiona,b,xreal.parameter:
m=0.7071!
m代表a=
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（3）=10.
b（ll）=20.
write（叮（13（13（f5・2』x"
）r）（（a（ij）沪1山3）启1J3）!
callagausfa^^xJJs）
wri9（^（3（5f8・2,/）r））（（:
输出结果endif
stopendsubroutineagaus（a,bnx」,js）
dimensiona（n,n）,x（n）,b（n）,js（n）
逻辑变fi：
d=abs（a（ij»
t=a（i,k）
a（i/k）=a（ijs（k））
a（i,js（k））=t
lf（is/=k）then
do
t=a（k,j）
a（k,j）=a（is,j）
writ屮JOO）
doj=k+l/n
a（i/j）=a（M）-a（i/kra（kJ）
b（i）=b（i）・a（i*广b（k）
write（*,100）
x（n）=b（n）/a（n/n）
do1=01,1,1
doj=i+l/n
100format（lx/faila）
dok=n/l/-l
endifenddoreturn
end精品文档word文档可以编借！
谢谢下载!
升级会员 