哈尔滨工程大学数值分析大作业附fortran程序.docx
- 文档编号:28585643
- 上传时间:2023-07-19
- 格式:DOCX
- 页数:22
- 大小:83.26KB
哈尔滨工程大学数值分析大作业附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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 哈尔滨工程 大学 数值 分析 作业 fortran 程序