北航数值分析报告大作业第三题fortran.docx

资源描述

北航数值分析报告大作业第三题fortran.docx

《北航数值分析报告大作业第三题fortran.docx》由会员分享，可在线阅读，更多相关《北航数值分析报告大作业第三题fortran.docx（32页珍藏版）》请在冰豆网上搜索。

北航数值分析报告大作业第三题fortran.docx

北航数值分析报告大作业第三题fortran

“数值分析“计算实习大作业第三题

——SY1415215

孔维鹏

一、计算说明

1、将

，

分别代入方程组（A.3）得到关于t，u，v，w的的方程组，调用离散牛顿迭代子函数求出与

，

对应的

，

。

2、调用分片二次代数插值子函数在点（

）处插值得到

），得到数表

。

3、对于

，分别调用最小二乘拟合子函数计算系数矩阵

及误差

，直到满足精度，即求得最小的k值及系数矩阵

。

4、将

，

分别代入方程组（A.3）得到关于

，

的的方程组，调用离散牛顿迭代子函数求出与

，

对应的

，

，调用分片二次代数插值子函数在点（

）处插值得到

）；调用步骤3中求得的系数矩阵

求得

，打印数表

。

二、源程序（FORTRAN）

PROGRAMSY1415215

DIMENSIONX（11）,Y（21）,T（6）,U（6）,Z（6,6）,UX（11,21）,TY（11,21）,FXY（11,21）,C（6,6）

DIMENSIONX1（8）,Y1（5）,FXY1（8,5）,PXY1（8,5）,UX1（8,5）,TY1（8,5）

REAL（8）X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1

OPEN（1,FILE='第三题计算结果.TXT'）

DOI=1,11

X（I）=0.08*（I-1）

ENDDO

DOI=1,21

Y（I）=0.5+0.05*（I-1）

ENDDO

*****求解非线性方程组，得到z=f（t，u）的函数*******

DOI=1,11

DOJ=1,21

CALLDISNEWTON_NONLINEAR（X（I）,Y（J）,UX（I,J）,TY（I,J））

ENDDO

*************分片二次插值得到z=f（x，y）***********

DOI=1,11

DOJ=1,21

CALLINTERPOLATION（UX（I,J）,TY（I,J）,FXY（I,J））

ENDDO

WRITE（1,"（'数表（x,y,f（x,y））:

'）"）

WRITE（1,"（3X,'X',7X,'Y',10X,'F（X,Y）'）"）

DOI=1,11

DOJ=1,21

WRITE（1,'（1X,F5.2,2X,F5.3,2X,E20.13）'）X（I）,Y（J）,FXY（I,J）

ENDDO

WRITE（1,"（''）"）

ENDDO

***********最小二乘拟合得到P（x，y）**************

N=11

M=21

WRITE（1,'（"","K和σ分别为：

"）'）

DOK=1,20

CALLLSFITTING（X,Y,FXY,C,N,M,K,K,E）

WRITE（1,'（I3,2X,E20.13）'）K-1,E

IF（E<10E-7）EXIT

ENDDO

WRITE（1,'（""）'）

WRITE（1,'（"系数矩阵Crs（按行）为：

"）'）

DOI=1,K

DOJ=1,K

WRITE（1,'（E20.13,2X,\）'）C（I,J）

ENDDO

WRITE（1,"（''）"）

WRITE（*,"（''）"）

ENDDO

DOI=1,8

X1（I）=0.1*I

ENDDO

DOJ=1,5

Y1（J）=0.5+0.2*J

ENDDO

DOI=1,8

DOJ=1,5

CALLDISNEWTON_NONLINEAR（X1（I）,Y1（J）,UX1（I,J）,TY1（I,J））

ENDDO

DOI=1,8

DOJ=1,5

CALLINTERPOLATION（UX1（I,J）,TY1（I,J）,FXY1（I,J））

ENDDO

PXY1=0

DOI=1,8

DOJ=1,5

DOII=1,K

DOJJ=1,K

PXY1（I,J）=PXY1（I,J）+C（II,JJ）*（X1（I）**（II-1））*（Y1（J）**（JJ-1））

ENDDO

WRITE（1,'（""）'）

WRITE（1,'（"数表（x,y,f（x,y）,p（x,y））:

"）'）

WRITE（1,"（2X,'X',6X,'Y',12X,'F（X,Y）',14X,'P（X,Y）'）"）

DOI=1,8

DOJ=1,5

WRITE（1,'（F5.3,2X,F5.3,2X,E20.13,2X,E20.13）'）X1（I）,Y1（J）,FXY1（I,J）,PXY1（I,J）

ENDDO

WRITE（1,"（''）"）

ENDDO

（1）

END

***********用离散牛顿法求解非线性方程组****************

SUBROUTINEDISNEWTON_NONLINEAR（X1,Y1,U,T）

PARAMETER（N=4）

REALEPS!

EPS为迭代精度，M为最大迭代次数

DIMENSIONX（N）,H（N）,Y（N）,JA（N,N）,E（N）,XK（N）

REAL（8）JA,X,H,Y,E,XK,U,T,V,W,X1,Y1,E1,E2

F1（T,U,V,W）=0.5*COS（T）+U+V+W-X1-2.67

F2（T,U,V,W）=T+0.5*SIN（U）+V+W-Y1-1.07

F3（T,U,V,W）=0.5*T+U+COS（V）+W-X1-3.74

F4（T,U,V,W）=T+0.5*U+V+SIN（W）-Y1-0.79

EPS=10E-12

M=100

X=1.0

DOK=1,M

H=1

计算Y=F（x）

（1）=F1（X

（1）,X

（2）,X（3）,X（4））

（2）=F2（X

（1）,X

（2）,X（3）,X（4））

Y（3）=F3（X

（1）,X

（2）,X（3）,X（4））

Y（4）=F4（X

（1）,X

（2）,X（3）,X（4））

计算JA（N,N）

E=0.0

DOI=1,N

DOJ=1,N

DOJJ=1,N

IF（JJ==J）THEN

E（JJ）=X（JJ）+H（JJ）

ELSE

E（JJ）=X（JJ）

ENDIF

ENDDO

IF（I==1）THEN

JA（I,J）=（F1（E

（1）,E

（2）,E（3）,E（4））-F1（X

（1）,X

（2）,X（3）,X（4）））/H（J）

ELSEIF（I==2）THEN

JA（I,J）=（F2（E

（1）,E

（2）,E（3）,E（4））-F2（X

（1）,X

（2）,X（3）,X（4）））/H（J）

ELSEIF（I==3）THEN

JA（I,J）=（F3（E

（1）,E

（2）,E（3）,E（4））-F3（X

（1）,X

（2）,X（3）,X（4）））/H（J）

ELSEIF（I==4）THEN

JA（I,J）=（F4（E

（1）,E

（2）,E（3）,E（4））-F4（X

（1）,X

（2）,X（3）,X（4）））/H（J）

ENDIF

ENDDO

求解线性方程组

CALLGAUSS（JA,XK,-Y,N）

判断精度

CALLNORM（XK,N,E1）

CALLNORM（X,N,E2）

IF（E1/E2<=EPS）THEN

T=X

（1）

U=X

（2）

EXIT

ELSE

DOI=1,N

X（I）=X（I）+XK（I）

ENDDO

ENDIF

ENDDO

RETURN

END

**********列主元高斯消去法求解线性方程组*********

SUBROUTINEGAUSS（A,X,B,N）

DIMENSIONA（N,N）,B（N）,X（N）,T（N,N）,TB（N）

REALM（N,N）

REAL（8）A,B,X,T

消元过程

DOK=1,N-1

TA=A（K,K）

TL=K

DOL=K+1,N

IF（（A（L,K）>TA）.OR.（A（L,K）==TA））THEN

TA=A（L,K）

TL=L

DOJ=K,N

T（K,J）=A（K,J）

A（K,J）=A（TL,J）

A（TL,J）=T（K,J）

ENDDO

TB（K）=B（K）

B（K）=B（TL）

B（TL）=TB（K）

ENDIF

ENDDO

DOI=K+1,N

M（I,K）=A（I,K）/A（K,K）

A（I,K）=0

DOJ=K+1,N

A（I,J）=A（I,J）-M（I,K）*A（K,J）

ENDDO

B（I）=B（I）-M（I,K）*B（K）

ENDDO

回代过程

X（N）=B（N）/A（N,N）

DOK=N-1,1,-1

S=0.0

DOJ=K+1,N

S=S+A（K,J）*X（J）

ENDDO

X（K）=（B（K）-S）/A（K,K）

ENDDO

RETURN

END

***********求向量的无穷数************

SUBROUTINENORM（X,N,A）

DIMENSIONX（N）

REAL（8）X,A

A=ABS（X

（1））

DOI=2,N

IF（ABS（X（I））>ABS（X（I-1）））THEN

A=ABS（X（I））

ENDIF

ENDDO

RETURN

END

**************分片二次代数插值**************

SUBROUTINEINTERPOLATION（U,V,W）

PARAMETER（N=6,M=6）

DIMENSIONX（N）,Y（M）,Z（M,N）,LK（3）,LR（3）

REAL（8）X,Y,Z,H,T

REAL（8）U,V,W,LK,LR!

U,V分别为插值点处的坐标，W为插值结果

INTEGERR

**********************数据赋值**********************

DATAY/0.0,0.2,0.4,0.6,0.8,1.0/

DATAX/0.0,0.4,0.8,1.2,1.6,2.0/

DATAZ/-0.5,-0.42,-0.18,0.22,0.78,1.5,&

&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&

&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&

&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&

&2.06,1.18,0.46,-0.1,-0.5,-0.74,&

&3.5,2.38,1.42,0.62,-0.02,-0.5/

H=0.4

T=0.2

******************计算K,R*************************

IF（U<=X

（2）+H/2）THEN

K=2

ELSEIF（U>X（N-1）-H/2）THEN

K=N-1

ELSE

DOI=3,N-2

IF（（U>X（I）-H/2）.AND.（U<=X（I）+H/2））THEN

K=I

ENDIF

ENDDO

ENDIF

IF（V<=Y

（2）+T/2）THEN

R=2

ELSEIF（V>Y（M-1）-T/2）THEN

R=M-1

ELSE

DOJ=3,M-2

IF（（V>Y（J）-T/2）.AND.（V<=Y（J）+T/2））THEN

R=J

ENDIF

ENDDO

ENDIF

I=K

J=R

（1）=（U-X（I））*（U-X（I+1））/（X（I-1）-X（I））/（X（I-1）-X（I+1））

（2）=（U-X（I-1））*（U-X（I+1））/（X（I）-X（I-1））/（X（I）-X（I+1））

LK（3）=（U-X（I））*（U-X（I-1））/（X（I+1）-X（I））/（X（I+1）-X（I-1））

（1）=（V-Y（J））*（V-Y（J+1））/（Y（J-1）-Y（J））/（Y（J-1）-Y（J+1））

（2）=（V-Y（J-1））*（V-Y（J+1））/（Y（J）-Y（J-1））/（Y（J）-Y（J+1））

LR（3）=（V-Y（J））*（V-Y（J-1））/（Y（J+1）-Y（J））/（Y（J+1）-Y（J-1））

W=0

DOK=1,3

DOR=1,3

W=W+LK（K）*LR（R）*Z（J+R-2,I+K-2）

ENDDO

RETURN

END

*******************最小二乘拟合子函数**************

SUBROUTINELSFITTING（X,Y,Z,A,N,M,P,Q,DT1）

INTEGERP,Q

DIMENSIONX（N）,Y（M）,Z（N,M）,A（P,Q）

DIMENSIONAPX（20）,APY（20）,BX（20）,BY（20）,U（20,20）,V（20,M）

DIMENSIONT（20）,T1（20）,T2（20）

REAL（8）X,Y,Z,A,DT1

DOI=1,P

DOJ=1,Q

A（I,J）=0.0

ENDDO

IF（P>N）P=N

IF（P>20）P=20

IF（Q>M）Q=M

IF（Q>20）Q=20

XX=0

YY=0

D1=N

APX

（1）=0.0

DOI=1,N

APX

（1）=APX

（1）+X（I）

ENDDO

APX

（1）=APX

（1）/D1

DOJ=1,M

V（1,J）=0.0

DOI=1,N

V（1,J）=V（1,J）+Z（I,J）

ENDDO

V（1,J）=V（1,J）/D1

ENDDO

IF（P>1）THEN

D2=0.0

APX

（2）=0.0

DOI=1,N

G=X（I）-APX

（1）

D2=D2+G*G

APX

（2）=APX

（2）+（X（I）-XX）*G*G

ENDDO

APX

（2）=APX

（2）/D2

（2）=D2/D1

DOJ=1,M

V（2,J）=0.0

DOI=1,N

G=X（I）-APX

（1）

V（2,J）=V（2,J）+Z（I,J）*G

ENDDO

V（2,J）=V（2,J）/D2

ENDDO

D1=D2

ENDIF

DOK=3,P

D2=0.0

APX（K）=0.0

DOJ=1,M

V（K,J）=0.0

ENDDO

DOI=1,N

G1=1.0

G2=X（I）-APX

（1）

DOJ=3,K

G=（X（I）-APX（J-1））*G2-BX（J-1）*G1

G1=G2

G2=G

ENDDO

D2=D2+G*G

APX（K）=APX（K）+X（I）*G*G

DOJ=1,M

V（K,J）=V（K,J）+Z（I,J）*G

ENDDO

DOJ=1,M

V（K,J）=V（K,J）/D2

ENDDO

APX（K）=APX（K）/D2

BX（K）=D2/D1

D1=D2

ENDDO

D1=M

APY

（1）=0.0

DOI=1,M

APY

（1）=APY

（1）+Y（I）

ENDDO

APY

（1）=APY

（1）/D1

DOJ=1,P

U（J,1）=0.0

DOI=1,M

U（J,1）=U（J,1）+V（J,I）

ENDDO

U（J,1）=U（J,1）/D1

ENDDO

IF（Q>1）THEN

D2=0.0

APY

（2）=0.0

DOI=1,M

G=Y（I）-APY

（1）

D2=D2+G*G

APY

（2）=APY

（2）+（Y（I））*G*G

ENDDO

APY

（2）=APY

（2）/D2

（2）=D2/D1

DOJ=1,P

U（J,2）=0.0

DOI=1,M

G=Y（I）-APY

（1）

U（J,2）=U（J,2）+V（J,I）*G

ENDDO

U（J,2）=U（J,2）/D2

ENDDO

D1=D2

ENDIF

DOK=3,Q

D2=0.0

APY（K）=0.0

DOJ=1,P

U（J,K）=0.0

ENDDO

DOI=1,M

G1=1.0

G2=Y（I）-APY

（1）

DOJ=3,K

G=（Y（I）-APY（J-1））*G2-BY（J-1）*G1

G1=G2

G2=G

ENDDO

D2=D2+G*G

APY（K）=APY（K）+Y（I）*G*G

DOJ=1,P

U（J,K）=U（J,K）+V（J,I）*G

ENDDO

DOJ=1,P

U（J,K）=U（J,K）/D2

ENDDO

APY（K）=APY（K）/D2

BY（K）=D2/D1

D1=D2

ENDDO

V（1,1）=1.0

V（2,1）=-APY

（1）

V（2,2）=1.0

DOI=1,P

DOJ=1,Q

A（I,J）=0.0

ENDDO

DOI=3,Q

V（I,I）=V（I-1,I-1）

V（I,I-1）=-APY（I-1）*V（I-1,I-1）+V（I-1,I-2）

IF（I>=4）THEN

DOK=I-2,2,-1

V（I,K）=-APY（I-1）*V（I-1,K）+V（I-1,K-1）-BY（I-1）*V（I-2,K）

ENDDO

ENDIF

V（I,1）=-APY（I-1）*V（I-1,1）-BY（I-1）*V（I-2,1）

ENDDO

DOI=1,P

IF（I==1）THEN

（1）=1.0

ELSEIF（I==2）THEN

（1）=-APX

（1）

（2）=1.0

（1）=T

（1）

（2）=T

（2）

ELSE

T（I）=T2（I-1）

T（I-1）=-APX（I-1）*T2（I-1）+T2（I-2）

IF（I>=4）THEN

DOK=I-2,2,-1

T（K）=-APX（I-1）*T2（K）+T2（K-1）-BX（I-1）*T1（K）

ENDDO

ENDIF

（1）=-APX（I-1）*T2

（1）-BX（I-1）*T1

（1）

T2（I）=T（I）

DOK=I-1,1,-1

T1（K）=T2（K）

T2（K）=T（K）

ENDDO

ENDIF

DOJ=1,Q

DOK=I,1,-1

DOL=J,1,-1

A（K,L）=A（K,L）+U（I,J）*T（K）*V（J,L）

ENDDO

DT1=0.0

DOI=1,N

X1=X（I）

DOJ=1,M

Y1=Y（J）

X2=1.0

DD=0.0

DOK=1,P

G=A（K,Q）

DOKK=Q-1,1,-1

G=G*Y1+A（K,KK）

ENDDO

G=G*X2

DD=DD+G

X2=X2*X1

ENDDO

DT=DD-Z（I,J）

DT1=DT1+DT*DT

ENDDO

RETURN

END

三、计算结果

数表（x,y,f（x,y））:

XYUXTYF（X,Y）

0.000.5001.3450.2430.17E+00

0.000.5501.3220.2690.66E+00

0.000.6001.2990.2950.35E+00

0.000.6501.2770.3220.94E+00

0.000.7001.2550.3500.30E-02

0.000.7501.2350.377-0.87E-01

0.000.8001.2150.406-0.58E+00

0.000.8501.1960.434-0.72E+00

0.000.9001.1770.463-0.54E+00

0.000.9501.1590.492-0.86E+00

0.001.0001.1420.521-0.01E+00

0.001.0501.1250.550-0.74E+00

0.001.1001.1090.580-0.06E+00

0.001.1501.0930.609-0.00E+00

0.001.2001.0790.639-0.18E+00

0.001.2501.0640.669-0.52E+00

0.001.3001.0500.699-0.19E+00

0.001.3501.0370.729-0.48E+00

0.001.4001.0240.759-0.68E+00

0.001.4501.0110.790-0.52E+00

0.001.5001.0000.820-0.29E+00

0.080.5001.4150.2280.67E+00

0.080.5501.3910.2530.08E+00

0.080.6001.3680.2790.02E+00

0.080.6501.3460.3060.47E+00

0.080.7001.3250.3330.57E+00

0.080.7501.3040.3600.48E-01

0.080.8001.2840.388-0.73E-01

0.080.8501.2650.416-0.16E+00

0.080.9001.2460.444-0.29E+00

0.080.9501.2290.473-0.36E+00

0.081.0001.2110.502-0.08E+00

0.081.0501.1940.531-0.29E+00

0.081.1001.1780.560-0.78E+00

0.081.1501.1630.589-0.93E+00

0.081.2001.1480.619-0.44E+00

0.081.2501.1330.649-0.92E+00

0.081.3001.1190.679-0.71E+00

0.081.3501.1060.709-0.98E+00

0.081.4001.0

展开阅读全文