中南大学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