哈尔滨工程大学数值分析大作业附fortran程序Word文件下载.docx
- 文档编号:16427907
- 上传时间:2022-11-23
- 格式:DOCX
- 页数:19
- 大小:83.12KB
哈尔滨工程大学数值分析大作业附fortran程序Word文件下载.docx
《哈尔滨工程大学数值分析大作业附fortran程序Word文件下载.docx》由会员分享,可在线阅读,更多相关《哈尔滨工程大学数值分析大作业附fortran程序Word文件下载.docx(19页珍藏版)》请在冰豆网上搜索。
二数学原理
对于方程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,'
)k,x(k)
if(abs(x(k)-x(k-1))<
1e-6)exit!
终止迭代条件
enddo
stop
end
(2)对于
,按照上述数学原理,编制的程序如下
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
)k,x(k)!
5e-6)exit!
四结果分析和讨论
(1)对于方程
,当初值分别为
时,所得结果如下
k
xk
初始值1
初始值2
初始值3
x0
1
0.45
0.65
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
-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
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
37
x37
结果分析与讨论:
从计算结果可以看到,当取的初始值不同时,虽然均得到了近似解x*=1.324718,但收敛速度明显不同。
当初始值x0充分接近于方程的单根时,可保证迭代序列快速收敛。
在本例中,初始值1、0.65和0.45距方程的单根越来越远,故收敛速度越来越慢。
(2)对于方程
,当初始值x0=2时计算结果如下
初始值
-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号铰接点有
programmain
integer,parameter:
n=13!
输入量:
系数阵的维数
js(n)
dimension:
a(n,n),b(n),x(n)
doubleprecisiona,b,x
real,parameter:
m=0.7071!
m代表α=
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
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
最大元素不在K行,K行
if(is/=k)then
t=a(k,j)
a(k,j)=a(is,j)
a(is,j)=t!
交换到K列
t=b(k)
b(k)=b(is)
b(is)=t
endif!
最大元素在主对角线上
消去
if(l==0)then
write(*,100)
return
doj=k+1,n
a(k,j)=a(k,j)/a(k,k)
b(k)=b(k)/a(k,k)!
求三角矩阵
doi=k+1,n
a(i,j)=a(i,j)-a(i,k)*a(k,j)
b(i)=b(i)-a(i,k)*b(k)
if(abs(a(n,n))+1.0==1.0)then
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)
x(i)=b(i)-t
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
return
end
由程序计算的各个杆的受力如下:
f1
f2
f3
f4
f5
f6
f7
-28.28
20.00
10.00
-30.00
14.14
0.00
f8
f9
f10
f11
f12
f13
7.07
25.00
-35.36
列方程时假定各杆均受拉力,得到的结果有负值,说明力与假设方向相反。
高斯消去法虽然能够执行,但随便用akk(k)作为主元,有时会扩大误差,导致结果不可靠,甚至严重失真,而高斯列主元消去法则不会有这种情况发生。
最初采用的是高斯-赛德尔迭代法解此矩阵,但是发现很难得到收敛解。
因为必须满足迭代条件才可以,而找到满足条件的迭代矩阵非常困难,故最终选取了没有收敛限制的直接法解此矩阵。
附录
题目一程序:
(1)
(2)
题目二程序:
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.,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./
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 哈尔滨工程 大学 数值 分析 作业 fortran 程序