中南大学fortran课程具体设计.docx
《中南大学fortran课程具体设计.docx》由会员分享，可在线阅读，更多相关《中南大学fortran课程具体设计.docx（19页珍藏版）》请在冰豆网上搜索。
中南大学fortran课程具体设计
文件
OPEN（10,FILE="X.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
DOI=1,10
X=10.0*I
WRITE（10,100,REC=I）X
100FORMAT（F5.1）
ENDDO
OPEN（20,FILE="Y.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
WRITE（20,200,REC=1）6.0
WRITE（20,200,REC=2）14.0
WRITE（20,200,REC=3）26.0
WRITE（20,200,REC=4）33.0
WRITE（20,200,REC=5）46.0
WRITE（20,200,REC=6）54.0
WRITE（20,200,REC=7）67.0
WRITE（20,200,REC=8）75.0
WRITE（20,200,REC=9）84.0
WRITE（20,200,REC=10）100.0
200FORMAT（F6.1）
END
DIMENSIONY（10）
OPEN（20,FILE="Y.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
DOI=1,10
READ（20,200,REC=I）Y（I）
ENDDO
X=10
H=10
EPS=1.0
DOT=15,95,10
CALLEEATK（X,H,10,Y,T,EPS,Z）
WRITE（*,10）T,Z
10FORMAT（1X,'压力P（KN）=',F5.1,10X,'变形（mm）=',F5.1）
200FORMAT（F6.1）
ENDDO
CLOSE（20）
END
SUBROUTINEEEATK（X1,H,N,Y,T,EPS,Z）
DIMENSIONY（N）,XM（10）,YM（10）
M=10
IF（M.GT.N）M=N
Z=0.0
IF（M.LE.0）RETURN
IF（N.EQ.1）THEN
Z=Y
（1）
RETURN
ENDIF
IF（M.EQ.1）M=2
IF（T.LE.X1）THEN
K=1
ELSEIF（T.GE.X1+（N-1）*H）THEN
K=N
ELSE
K=1
J=N
10IF（IABS（K-J）.NE.1）THEN
L=（K+J）/2
IF（T.LT.X1+（L-1）*H）THEN
J=L
ELSE
K=L
ENDIF
GOTO10
ENDIF
IF（ABS（T-X1-（L-1）*H）.GT.ABS（T-X1-（J-1）*H））K=J
ENDIF
J=1
L=0
DO20I=1,M
K=K+J*L
IF（（K.LT.1）.OR.（K.GT.N））THEN
L=L+1
J=-J
K=K+J*L
ENDIF
XM（I）=X1+（K-1）*H
YM（I）=Y（K）
L=L+1
J=-J
20CONTINUE
I=2
P=1.0+EPS
100IF（（I.LE.M）.AND.（P.GE.EPS））THEN
Z=YM（I）
DO30J=2,I
30Z=YM（J-1）+（T-XM（J-1））*（YM（J-1）-Z）/（XM（J-1）-XM（I））
YM（I）=Z
P=ABS（YM（I）-YM（I-1））
I=I+1
GOTO100
ENDIF
RETURN
END
切比雪夫曲线拟合
DIMENSIONX（10）,Y（10）,A（7）
N=10
M=6
OPEN（10,FILE="X.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
DOI=1,10
READ（10,100,REC=I）X（I）
OPEN（20,FILE="Y.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
READ（20,200,REC=I）Y（I）
ENDDO
10CONTINUE
CALLHCHIR（X,Y,N,A,M,M1）
WRITE（*,20）（I,A（I）,I=1,M）
WRITE（*,*）
20FORMAT（1X,'A（',I2,'）=',F15.6）
100FORMAT（F5.1）
200FORMAT（F6.1）
OPEN（30,FILE="A.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=20）
DOI=1,6
WRITE（30,300,REC=I）A（I）
ENDDO
300FORMAT（F15.6）
END
SUBROUTINEHCHIR（X,Y,N,A,M,M1）
DIMENSIONX（N）,Y（N）,A（M1）,IX（20）,H（20）
DO5I=1,M1
5A（I）=0.0
IF（M.GE.N）M=N-1
IF（M.GE.20）M=19
M1=M+1
HA=0.0
IX
（1）=1
IX（M1）=N
L=（N-1）/M
J=L
DO10I=2,M
IX（I）=J+1
J=J+L
10CONTINUE
20HH=1.0
DO30I=1,M1
A（I）=Y（IX（I））
H（I）=-HH
HH=-HH
30CONTINUE
DO50J=1,M
II=M1
Y2=A（II）
H2=H（II）
DO40I=J,M
D=X（IX（II））-X（IX（M1-I））
Y1=A（M-I+J）
H1=H（M-I+J）
A（II）=（Y2-Y1）/D
H（II）=（H2-H1）/D
II=M-I+J
Y2=Y1
H2=H1
40CONTINUE
50CONTINUE
HH=-A（M1）/H（M1）
DO60I=1,M1
60A（I）=A（I）+H（I）*HH
DO80J=1,M-1
II=M-J
D=X（IX（II））
Y2=A（II）
DO70K=M1-J,M
Y1=A（K）
A（II）=Y2-D*Y1
Y2=Y1
II=K
70CONTINUE
80CONTINUE
HM=ABS（HH）
IF（HM.LE.HA）THEN
A（M1）=-HM
RETURN
ENDIF
A（M1）=HM
HA=HM
IM=IX
（1）
H1=HH
J=1
DO100I=1,N
IF（I.EQ.IX（J））THEN
IF（J.LT.M1）J=J+1
ELSE
H2=A（M）
DO90K=M-1,1,-1
90H2=H2*X（I）+A（K）
H2=H2-Y（I）
IF（ABS（H2）.GT.HM）THEN
HM=ABS（H2）
H1=H2
IM=I
ENDIF
ENDIF
100CONTINUE
IF（IM.EQ.IX
（1））RETURN
I=1
110IF（IM.GE.IX（I））THEN
I=I+1
IF（I.LE.M1）GOTO110
ENDIF
IF（I.GT.M1）I=M1
IF（I.EQ.（I/2）*2）THEN
H2=HH
ELSE
H2=-HH
ENDIF
IF（H1*H2.GE.0.0）THEN
IX（I）=IM
GOTO20
ENDIF
IF（IM.LT.IX
（1））THEN
DO120J=M,1,-1
120IX（J+1）=IX（J）
IX
（1）=IM
GOTO20
ENDIF
IF（IM.GT.IX（M1））THEN
DO130J=2,M1
130IX（J-1）=IX（J）
IX（M1）=IM
GOTO20
ENDIF
IX（I-1）=IM
GOTO20
END
PROGRAMMAIN
DIMENSIONX（10）,Y（10）,S（9）
REAL（8）,ALLOCATABLE:
:
S1（:
:
）,S2（:
）,A（:
）,JS（:
）
REAL（8）X,Y,S,R
N=6
ALLOCATE（S1（N+1,N+1）,S2（N+1）,A（N+1）,JS（N+1））
K=1
OPEN（10,FILE="X.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
DOI=1,10
READ（10,300,REC=I）X（I）
K=K+1
ENDDO
CLOSE（10）
K=1
OPEN（20,FILE="Y.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=10）
READ（20,400,REC=I）Y（I）
K=K+1
CLOSE（20）
300FORMAT（F5.1）
400FORMAT（F6.1）
I=1
DOI=1,9
S（I）=5+10*I
ENDDO
WRITE（*,200）'压力（KN）','变形（mm）'
DOI=1,9
CALLLAGRANGE（X,Y,10,S（I）,R）
WRITE（*,100）S（I）,R
ENDDO
200FORMAT（1X,A8,10X,A8）
100FORMAT（1x,F8.3,10X,F8.3）
DOI=1,N+1
DOJ=1,N+1
CALLSUM1（S1（I,J）,X,I+J-2）
ENDDO
ENDDO
S1（1,1）=K+1
DOI=1,N+1
CALLSUM2（S2（I）,X,Y,I-1）
ENDDO
CALLGUASS（S1,S2,N+1,A,JS）
OPEN（30,FILE="A.TXT",FORM="FORMATTED",ACCESS="DIRECT",RECL=20）
DOI=1,6
WRITE（30,500,REC=I）A（I）
ENDDO
500FORMAT（F15.6）
WRITE（*,*）'拟合曲线的系数为:
',A
END
SUBROUTINESUM1（S,X,N）
DIMENSIONX（10）
REAL（8）X,S
DOI=1,10
S=S+X（I）**N
ENDDO
END
SUBROUTINELAGRANGE（X,Y,N,S,R）
DIMENSIONX（N）,Y（N）
REAL（8）:
:
X,Y,S,R
R=0.0
IF（N<=0）THEN
PRINT*,"ERROR"
ENDIF
IF（N==1）R=Y
（1）
IF（N==2）R=（S-X
（2））/（Y
（1）-Y
（2））+（S-X
（1））/（Y
（2）-Y
（1））
I=1
DOWHILE（x（i）<=s）
I=I+1
ENDDO
K=I-4
if（K<1）K=1
M=I+3
IF（M>=N）M=N
DOJ=K,M
T=1.0
DOL=K,M
IF（L/=J）T=T*（S-X（L））/（X（J）-X（L））
ENDDO
R=R+T*Y（J）
ENDDO
END
SUBROUTINESUM2（S,X,Y,N）
DIMENSIONX（10）,Y（10）
REAL（8）X,Y,S
DOI=1,10
S=S+Y（I）*X（I）**N
ENDDO
END
subroutineGuass（a,b,n,x,js）
real（8）a（n,n）,b（n）,x（n）,js（n）
real（8）T,D
!
全选主元法，从右下角n-k+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==0.0）then
write（*,*）'fail'
stop
else
!
行调整
if（is/=k）then
doj=k,n
T=a（k,j）
a（k,j）=a（is,j）
a（is,j）=T
enddo
!
对b数组做出调整
T=b（k）
b（k）=b（is）
b（is）=T
endif
!
对行进行变换
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
endif
endif
!
高斯消元法的核心变换
doi=k+1,n
a（i,k）=a（i,k）/a（k,k）
enddo
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!
结束最外层循环
!
回代
x（n）=b（n）/a（n,n）
doi=n-1,1,-1
T=0
doj=i+1,n
T=T+a（i,j）*x（j）
enddo
x（i）=（b（i）-T）/a（i,i）
enddo
!
由于在全主元法中进行了列变换，所以最后应对x（i）进行对应的调整
dok=n-1,1,-1
if（js（k）/=k）then
T=x（k）
x（k）=x（js（k））
x（js（k））=T
endif
enddo
end
图像
usemsflib
integerstatus
type（xycoord）xy
draw=setcolorrgb（#00ff00）
callmoveto（40,150,xy）
calloutgtext（"0"）
callmoveto（25,50,xy）
calloutgtext（"100"）
callmoveto（30,100,xy）!
画标注
calloutgtext（"50"）
callmoveto（160,150,xy）
calloutgtext（"变形/mm"）
callmoveto（120,20,xy）
calloutgtext（"绘制压力-变形拟合曲线"）
x=50
y=50
callmoveto（50,50,xy）
wx=-5
wy=5
x=x+wx
y=y+wy
draw=lineto（x,y）
callmoveto（55,55,xy）!
纵坐标画箭头
wx=5
wy=-5
x=x+wx
y=y+wy
draw=lineto（x,y）
status=setwindow（.true.,-55,-50,30,20）
callmoveto（50,50,xy）
wx=0
wy=100
x=x+wx
y=y+wy
draw=lineto（x,y）!
画坐标轴
wx=100
wy=0
x=x+wx
y=y+wy
draw=lineto（x,y）
wx=-5
wy=-5
x=x+wx
y=y+wy
draw=lineto（x,y）!
横坐标轴画箭头
callmoveto（145,155,xy）
wx=5
wy=5
x=x+wx
y=y+wy
draw=lineto（x,y）
dox=0,100,0.0001
y=-11.639880+2.106925*X-0.060745*X**2+0.001417**3-0.000015**4
status=setpixelrgb_w（x,y,#ffff00）
enddo
end
USEMSFLIB
TYPE（WXYCOORD）:
:
T
INTEGER:
:
RESULT
REAL（8）X,Y
CHARACTER*4STR
RESULT=SETWINDOW（.TRUE.,-40.,-40.,150.,150.）
RESULT=SETCOLORRGB（#00FF00）
CALLMOVETO_W（-1.05D1,0D0,T）
RESULT=LINETO_W（1.05D2,0D0）
CALLMOVETO_W（0D0,-1.05D1,T）
RESULT=LINETO_W（0D0,1.05D2）
calloutgtext（"绘制压力-变形拟合曲线"）
DOX=0,100,0.001
y=x
RESULT=SETPIXEL_W（X,Y）
ENDDO
DOX=0,100,10
CALLMOVETO_W（X-0.3,0D0,T）
WRITE（STR,"（F4.0）"）X
CALLOUTGTEXT（STR（1:
3））
ENDDO
DOy=0,100,10
CALLMOVETO_W（0D0,y-0.3,T）
WRITE（STR,"（F4.0）"）y
CALLOUTGTEXT（STR（1:
3））
ENDDO
END
升级会员 