终止迭代条件
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=njif(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文档可以编借!
谢谢下载!