数值计算基础上机答案.docx
《数值计算基础上机答案.docx》由会员分享，可在线阅读，更多相关《数值计算基础上机答案.docx（16页珍藏版）》请在冰豆网上搜索。
数值计算基础上机答案
Exe1
programmain
implicitnone
integer:
:
i
real（8）:
:
x（5）,y（5）,s（5）,t（5）,Lagrange,Newton,xx
data（x（i）,i=1,5）/1.615,1.634,1.702,1.828,1.921/
data（y（i）,i=1,5）/2.4145,2.46459,2.65271,3.03035,3.34066/
xx=1.682
print*,Lagrange（2,x（1:
3）,y（1:
3）,xx）
print*,Newton（2,x（1:
3）,y（1:
3）,xx）
data（s（i）,i=1,5）/0,0.1,0.2,0.3,0.4/
data（t（i）,i=1,5）/1,0.995,0.98007,0.95534,0.92106/
xx=0.048
print*,Lagrange（4,s,t,xx）
print*,Newton（4,s,t,xx）
endprogrammain
!
Lagrange插值
!
参数：
n插值多项式阶数，（x,y）插值节点，xx插值点
functionLagrange（n,x,y,xx）
implicitnone
integer:
:
n,i,j
real（8）,intent（in）:
:
x（0:
n）,y（0:
n）,xx
real（8）:
:
Lagrange,l
Lagrange=0
doi=0,n
l=y（i）
doj=0,n
if（j.ne.i）then
l=l*（xx-x（j））/（x（i）-x（j））
endif
enddo
Lagrange=Lagrange+l
enddo
endfunctionLagrange
!
Newton插值
!
参数：
n插值多项式阶数，（x,y）插值节点，xx插值点
functionNewton（n,x,y,xx）
implicitnone
integer:
:
n,i,j
real（8）,intent（in）:
:
x（0:
n）,y（0:
n）,xx
real（8）:
:
Newton,P（0:
n,0:
n）,Q
P（:
0）=y
Newton=P（0,0）
doi=1,n
doj=i,n
P（j,i）=（P（j,i-1）-P（j-1,i-1））/（x（j）-x（j-i））
enddo
Q=1
doj=1,i
Q=Q*（xx-x（j-1））
enddo
Newton=Newton+Q*P（i,i）
enddo
endfunctionNewton
Exe2
programmain
implicitnone
real,EXTERNAL:
:
f,g
realT1,T2,Romberg,R
print*,T1（f,0.0,1.0,8）
print*,T2（f,0.0,1.0,8）
print*,Romberg（f,0.0,1.0,3）
print*,Romberg（g,0.0,1.0,4）
print*,R（f,0.0,1.0,3,3）
print*,R（g,0.0,1.0,4,3）
endprogrammain
realfunctionf（x）
implicitnone
realx
if（x.eq.0）then
f=1
else
f=sin（x）/x
endif
endfunctionf
REALFUNCTIONg（x）
IMPLICITNONE
realx
g=4/（1+x**2）
endFUNCTIONg
!
复合梯形求积
realfunctionT1（f,a,b,n）
implicitnone
real,EXTERNAL:
:
f
realx,a,b
integeri,n
T1=f（a）+f（b）
doi=1,n-1
x=a+（b-a）*i/n
T1=T1+2*f（x）
enddo
T1=T1*（b-a）/n/2
endfunctionT1
!
梯形求积逐次对分递推
recursiverealfunctionT2（f,a,b,k）result（T2result）
implicitnone
real,EXTERNAL:
:
f
reala,b,y
integeri,j,k
if（k.eq.0）then
T2result=（f（a）+f（b））*（b-a）/2
else
y=0
doj=1,2**（k-1）
y=y+f（a+（2*j-1）*（b-a）/2**k）
enddo
T2result=T2（f,a,b,k-1）/2+（b-a）*y/2**k
endif
endfunctionT2
!
Simpson求积
REALFUNCTIONSimpson（f,a,b,k）
IMPLICITNONE
real,EXTERNAL:
:
f
reala,b,T2
INTEGERk
if（k.ge.1）THEN
Simpson=（4*T2（f,a,b,k）-T2（f,a,b,k-1））/3
else
PRINT*,"Error:
degreeofSimpsonintegralmustbe>=1"
STOP
endif
endFUNCTIONSimpson
!
Cotes求积
REALFUNCTIONCotes（f,a,b,k）
IMPLICITNONE
real,EXTERNAL:
:
f
reala,b,Simpson
INTEGERk
if（k.ge.2）THEN
Cotes=（16*Simpson（f,a,b,k）-Simpson（f,a,b,k-1））/15
else
PRINT*,"Error:
degreeofSimpsonintegralmustbe>=2"
STOP
endif
endFUNCTIONCotes
!
Romberg求积
REALFUNCTIONRomberg（f,a,b,k）
IMPLICITNONE
real,EXTERNAL:
:
f
reala,b,Cotes
INTEGERk
if（k.ge.3）THEN
Romberg=（64*Cotes（f,a,b,k）-Cotes（f,a,b,k-1））/63
else
PRINT*,"Error:
degreeofSimpsonintegralmustbe>=3"
STOP
endif
endFUNCTIONRomberg
!
Romberg求积递推
!
参数列表：
f被积函数,a积分下限,b积分上限,k递推次数,m加速次数
recursiveREALFUNCTIONR（f,a,b,k,m）result（RR）
IMPLICITNONE
real,EXTERNAL:
:
f
reala,b,y
INTEGERk,m,j
if（k.ge.m）THEN
if（m.gt.0）THEN
RR=（4**m*R（f,a,b,k,m-1）-R（f,a,b,k-1,m-1））/（4**m-1）
ELSEIF（m.eq.0）THEN
if（k.eq.0）then
RR=（f（a）+f（b））*（b-a）/2
else
y=0
doj=1,2**（k-1）
y=y+f（a+（2*j-1）*（b-a）/2**k）
enddo
RR=R（f,a,b,k-1,m）/2+（b-a）*y/2**k
endif
else
PRINT*,"Error:
orderofRombergintegralmustbe>=0"
STOP
endif
else
PRINT*,"Error:
Rombergintegralmustmatchdegree>=order"
STOP
endif
endFUNCTIONR
Exe3
PROGRAMMAIN
IMPLICITNONE
real,EXTERNAL:
:
f
realDichotomy,Secant
PRINT*,Dichotomy（f,1e-6,1.0,2.0）
PRINT*,Secant（f,1e-6,1.0,2.0）
ENDPROGRAMMAIN
!
待解方程
REALFUNCTIONf（x）
realx
f=x**3-x-1
ENDFUNCTIONf
!
二分法
realFUNCTIONDichotomy（f,limit,a,b）
implicitnone
real,EXTERNAL:
:
f
reala,b,limit,xa,xb
IF（SIGN（1.0,f（a））*SIGN（1.0,f（b））.gt.0）THEN
PRINT*,"Error:
f（a）*f（b）isgreaterthanzero."
ENDIF
xa=a
xb=b
dowhile（abs（xa-xb）.ge.limit）
Dichotomy=（xa+xb）/2
IF（SIGN（1.0,f（Dichotomy））*SIGN（1.0,f（xa））.lt.0）THEN
xb=Dichotomy
ELSE
xa=Dichotomy
ENDIF
enddo
Dichotomy=（xa+xb）/2
ENDFUNCTIONDichotomy
!
弦截法
REALFUNCTIONSecant（f,limit,a,b）
IMPLICITNONE
real,EXTERNAL:
:
f
reallimit,a,b,e,x0,x1,x2
x0=a
x1=b
e=limit
doWHILE（e.ge.limit）
x2=x1-f（x1）*（x1-x0）/（f（x1）-f（x0））
e=abs（x2-x1）
x0=x1
x1=x2
enddo
Secant=x2
endFUNCTIONSecant
Exe4
programMain
implicitnone
REAL（8）:
:
A1（3,3）,B1（3）,A2（3,3）,B2（3）,A3（2,2）,B3
（2）,limit
integer:
:
i,j
REAL（8）:
:
x（3）,y
（2）
limit=1e-12
data（（A1（i,j）,j=1,3）,i=1,3）/7,1,2,2,8,2,2,2,9/
data（B1（i）,i=1,3）/10,8,6/
CALLGauss（3,A1,B1,x,.False.）
PRINT"（3f18.15）",x
CALLGauss（3,A1,B1,x,.TRUE.）
PRINT"（3f18.15）",x
x=（/0,0,0/）
CALLIteration（3,A1,B1,x,limit,0）
PRINT"（3f18.15）",x
x=（/0,0,0/）
CALLIteration（3,A1,B1,x,limit,1）
PRINT"（3f18.15）",x
CALLDoolittle（3,A1,b1,x）
PRINT"（3f18.15）",x
data（（A2（i,j）,j=1,3）,i=1,3）/4,0.24,-0.08,0.09,3,-0.15,0.04,-0.08,4/
data（B2（i）,i=1,3）/8,9,10/
CALLGauss（3,A2,B2,x,.False.）
PRINT"（3f18.15）",x
CALLGauss（3,A2,B2,x,.TRUE.）
PRINT"（3f18.15）",x
x=（/0,0,0/）
CALLIteration（3,A2,B2,x,limit,0）
PRINT"（3f18.15）",x
x=（/0,0,0/）
CALLIteration（3,A2,B2,x,limit,1）
PRINT"（3f18.15）",x
CALLDoolittle（3,A2,b2,x）
PRINT"（3f18.15）",x
!
data（（A3（i,j）,j=1,2）,i=1,2）/1e-9,1,1,1/
!
data（B3（i）,i=1,2）/1,2/
!
CALLGauss（2,A3,B3,y,.False.）
!
PRINT*,y
!
CALLGauss（2,A3,B3,y,.TRUE.）
!
PRINT*,y
!
y=（/0,0/）
!
CALLIteration（2,A3,B3,y,limit,0）
!
PRINT*,y
!
y=（/0,0/）
!
CALLIteration（2,A3,B3,y,limit,1）
!
PRINT*,y
!
CALLDoolittle（2,A3,b3,y）
!
PRINT*,y
endprogramMain
!
Gauss消去法
!
参数：
n方程阶数，A系数矩阵，b右侧向量，x解向量，opt是否列主元
SUBROUTINEGauss（n,A,b,x,opt）
IMPLICITNONE
integer:
:
n,i,j,k
REAL（8）:
:
x（n）
REAL（8）,intent（in）:
:
A（n,n）,b（n）
REAL（8）:
:
Ab（n,n+1）,swap（n+1）,temp
LOGICAL,intent（in）:
:
opt
!
增广矩阵
Ab（1:
n,1:
n）=A
Ab（1:
n,n+1）=b
doj=1,n-1
!
doi=1,n
!
WRITE（*,"（1x,4f9.2）"）Ab（i,:
）
!
enddo
!
列主元
if（opt）then
k=MAXLOC（abs（Ab（j:
n,j））,1）+j-1
if（k.gt.j）THEN
swap=Ab（k,:
）
Ab（k,:
）=Ab（j,:
）
Ab（j,:
）=swap
endif
endif
!
doi=1,n
!
WRITE（*,"（1x,4f9.2）"）Ab（i,:
）
!
enddo
!
print*,"-------------"
!
消去
doi=j+1,n
temp=Ab（i,j）/Ab（j,j）
Ab（i,:
）=Ab（i,:
）-temp*Ab（j,:
）
enddo
enddo
!
回带
x（n）=Ab（n,n+1）/Ab（n,n）
doi=n-1,1,-1
x（i）=Ab（i,n+1）
doj=i+1,n
x（i）=x（i）-Ab（i,j）*x（j）
enddo
x（i）=x（i）/Ab（i,i）
enddo
endSUBROUTINEGauss
!
迭代
!
参数：
n方程阶数，A系数矩阵，b右侧向量，x解向量（输入初始解，输出迭代解），limit误差限，opt选项（0：
Jacobi，1：
Gauss-Seidel）
SUBROUTINEIteration（n,A,b,x,limit,opt）
IMPLICITNONE
integer:
:
n,i,j
integer,INTENT（in）:
:
opt
REAL（8）:
:
x（n）,xx（n）,xSeidel（n）,e
REAL（8）,intent（in）:
:
A（n,n）,b（n）,limit
e=limit
dowhile（e.ge.limit）
doi=1,n
!
Gauss-Seidel迭代，保存现场
if（opt.eq.1）THEN
xSeidel=x
endif
!
迭代
xx（i）=b（i）
doj=1,n
if（i.ne.j）THEN
xx（i）=xx（i）-A（i,j）*x（j）
endif
enddo
xx（i）=xx（i）/A（i,i）
!
Gauss-Seidel迭代，代入新值
if（opt.eq.1）THEN
x（i）=xx（i）
endif
!
范数、终止依据
if（i.eq.n）THEN
SELECTCASE（opt）
CASE（0）
e=MAXVAL（abs（x-xx））
x=xx
!
WRITE（*,"（1x,3f12.9）"）x
CASE
（1）
e=MAXVAL（abs（xSeidel-x））
!
WRITE（*,"（1x,3f12.9）"）x
endSELECT
endif
enddo
enddo
ENDSUBROUTINEIteration
!
Doolittle分解
!
参数：
n方程阶数，A系数矩阵，b右侧向量，x解向量（输入初始解，输出迭代解）
SUBROUTINEDoolittle（n,A,b,x）
IMPLICITNONE
integer:
:
n,i,j
REAL（8）:
:
x（n）,L（n,n）,U（n,n）,y（n）
REAL（8）,intent（in）:
:
A（n,n）,b（n）
!
LU分解
U=0
L=0
U（1,:
）=A（1,:
）
L（:
1）=A（:
1）/U（1,1）
y
（1）=b
（1）
doi=2,n
U（i,i）=A（i,i）-SUM（L（i,1:
i-1）*U（1:
i,i））
L（i,i）=1
doj=i+1,n
U（i,j）=A（i,j）-SUM（L（i,1:
i-1）*U（1:
i,j））
L（j,i）=（A（j,i）-SUM（L（j,1:
i-1）*U（1:
i-1,i）））/U（i,i）
enddo
y（i）=b（i）-SUM（L（i,1:
i-1）*y（1:
i-1））
enddo
!
求解x
x（n）=y（n）/U（n,n）
doi=n-1,1,-1
x（i）=（y（i）-SUM（U（i,i+1:
n）*x（i+1:
n）））/U（i,i）
enddo
!
显示中间变量L、U、y
!
print*,"L"
!
doi=1,n
!
WRITE（*,"（1x,3f9.5）"）L（i,:
）
!
enddo
!
print*,"U"
!
doi=1,n
!
WRITE（*,"（1x,3f9.5）"）U（i,:
）
!
enddo
!
print*,"y"
!
print*,y
endSUBROUTINEDoolittle
升级会员 