哈尔滨工程大学数值分析资料报告大作业附fortran程序文档格式.docx
- 文档编号:16691166
- 上传时间:2022-11-25
- 格式:DOCX
- 页数:18
- 大小:81.61KB
哈尔滨工程大学数值分析资料报告大作业附fortran程序文档格式.docx
《哈尔滨工程大学数值分析资料报告大作业附fortran程序文档格式.docx》由会员分享,可在线阅读,更多相关《哈尔滨工程大学数值分析资料报告大作业附fortran程序文档格式.docx(18页珍藏版)》请在冰豆网上搜索。
。
当选择初值
时给出结果并分析现象,当
,迭代停止。
二数学原理
对于方程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〕对于
,按照上述数学原理,编制的程序如下
初始值x(0)
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(*,'
)k,x(k)!
write(10,'
if(abs(x(k)-x(k-1))<
5e-6)exit!
四结果分析和讨论
〔1〕对于方程
,当初值分别为
时,所得结果如下
k
xk
初始值1
初始值2
初始值3
x0
1
x1
2
x2
3
x3
4
x4
5
x5
6
x6
7
x7
8
x8
9
x9
10
x10
11
x11
12
x12
13
x13
14
x14
15
x15
16
x16
17
x17
18
x18
19
x19
20
x20
21
x21
22
x22
23
x23
24
x24
25
x25
26
x26
27
x27
28
x28
29
x29
30
x30
31
x31
32
x32
33
x33
34
x34
35
x35
36
x36
37
x37
结果分析与讨论:
从计算结果可以看到,当取的初始值不同时,虽然均得到了近似解x*=1.324718,但收敛速度明显不同。
当初始值x0充分接近于方程的单根时,可保证迭代序列快速收敛。
在本例中,初始值1、0.65和0.45距方程的单根越来越远,故收敛速度越来越慢。
〔2〕对于方程
,当初始值x0=2时计算结果如下
初始值
牛顿法有明显的几何解释,方程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〕做初等行变换使原方程组转化为:
按x4→x3→x2→x1的顺序回代求解出方程组的解。
针对各个铰接点列方程:
对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
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
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
f8
f9
f10
f11
f12
f13
列方程时假定各杆均受拉力,得到的结果有负值,说明力与假设方向相反。
高斯消去法虽然能够执行,但随便用akk(k)作为主元,有时会扩大误差,导致结果不可靠,甚至严重失真,而高斯列主元消去法如此不会有这种情况发生。
最初采用的是高斯-赛德尔迭代法解此矩阵,但是发现很难得到收敛解。
因为必须满足迭代条件才可以,而找到满足条件的迭代矩阵非常困难,故最终选取了没有收敛限制的直接法解此矩阵。
附录
题目一程序:
〔1〕
〔2〕
题目二程序:
m代表α=
方程的系数阵
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.,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 程序