?? pre_statically_load.f90
字號:
W=DOT_PRODUCT(PFD,UW)
WRITE(8,50)U(1),W
IF(U(1)<T) GOTO 60
50 FORMAT(1X,E15.5,1X,E15.8)
100 FORMAT(11(1X,E15.8))
CONTAINS
!最速下降法解非線性方程組子程序
SUBROUTINE NLGRAD(NN,X,C,CX,DF,EPS,MM)
DOUBLE PRECISION X(NN),C,CX(NN),DF(NN),EPS,ASF,ASFH
3 DO I=1,NN
IF(X(I)==0.0)THEN
CX(I)=C
ELSE
CX(I)=C*X(I)
ENDIF
ENDDO
ASF=FNC(NN,X)
WRITE(*,*)ASF
IF(ASF<EPS) RETURN
SUM=0.0
DO I=1,NN
X(I)=X(I)+CX(I)
ASFH=FNC(NN,X)
DF(I)=(ASFH-ASF)/CX(I)
SUM=SUM+DF(I)**2
X(I)=X(I)-CX(I)
ENDDO
RLTA=ASF/SUM
DO I=1,NN
X(I)=X(I)-RLTA*DF(I)
ENDDO
MM=MM+1
GOTO 3
END SUBROUTINE NLGRAD
!預靜載非線性方程組左端函數值子程序
FUNCTION FNC(NN,ASX)
DOUBLE PRECISION ASX(NN),ASY(NN),ASS(N*M,1),FNC
FNC=0.0
DO J=1,N*M
ASS(J,1)=ASX(J)
ENDDO
DO J=1,N*M
LLL2(J,1)=ASS(J,1)-A0(J,1)
ENDDO
CALL INTEGRAL2(GG1,P10,F10,ASS,P10,F10,ASS,GM0,CI0,N,M,HX,HY)
CALL INTEGRAL2(GG2,P0,F20,ASS,P20,F0,ASS,GM0,CI0,N,M,HX,HY)
CALL INTEGRAL2(GG6,P10,F10,A0,P10,F10,A0,GM0,CI0,N,M,HX,HY)
CALL INTEGRAL2(GG7,P0,F20,A0,P20,F0,A0,GM0,CI0,N,M,HX,HY)
DO J=1,N*M
LL1(J,1)=E*(GG1(J,1)-GG2(J,1)-GG6(J,1)+GG7(J,1))
ENDDO
FS=MATMUL(M2N,LL1)
CALL INTEGRAL2(GG3,GM20,CI0,FS,P0,F20,ASS,P0,F0,N,M,HX,HY)
CALL INTEGRAL2(GG4,GM0,CI20,FS,P20,F0,ASS,P0,F0,N,M,HX,HY)
CALL INTEGRAL2(GG5,GM10,CI10,FS,P10,F10,ASS,P0,F0,N,M,HX,HY)
CALL INTEGRAL4(GG8,GMB2,CI0,FS,PB,F20,ASS,PB,F0,N,M,MK,HX)
CALL INTEGRAL4(GG9,GMB,CI20,FS,PB,F20,ASS,PB,F0,N,M,MK,HX)
CALL INTEGRAL4(GG10,PB,F20,LLL2,PB,F20,ASS,PB,F0,N,M,MK,HX)
LL2=MATMUL(M3,LLL2)
LL3=MATMUL(SS5,FS)
LL4=MATMUL(SS6,FS)
LL5=MATMUL(SS3,LLL2)
LL6=MATMUL(SS2,LLL2)
LL7=MATMUL(S5,ASS)
LL8=MATMUL(SS4,ASS)
DO J=1,N*M
ASY(J)=LL2(J,1)&
+PP0*LL7(J,1)&
-HP*(GG3(J,1)+GG4(J,1)-2.0*GG5(J,1))&
-(EK*AK*EEK/E*(LL3(J,1)-VK*LL4(J,1))&
-EK*IXK*LL5(J,1)&
-GK*JK*LL6(J,1)&
-EK*AK/E/HP*PP0*LL8(J,1)&
+EK*AK/E*(GG8(J,1)-VK*GG9(J,1))&
-EK*AK*EEK*GG10(J,1))&
-Q0*QQ0(J,1)
ENDDO
DO J=1,N*M
FNC=FNC+ASY(J)**2
ENDDO
RETURN
END FUNCTION FNC
!Runge-Kutta法子程序
SUBROUTINE RGKT2(NN,DT,U,DU,UC,UU,N,M,MK)
DOUBLE PRECISION U(NN),DU(NN),UC(NN),UU(NN),SN(4),DT
SN(1)=0.5*DT
SN(2)=SN(1)
SN(3)=DT
SN(4)=DT
DO I=1,NN
UU(I)=U(I)
ENDDO
DO J1=1,3
DO I=1,NN
UC(I)=UU(I)+SN(J1)*DU(I)
U(I)=U(I)+SN(J1+1)*DU(I)/3.0
ENDDO
CALL DIF(NN,UC,DU,N,M,MK)
ENDDO
DO I=1,NN
U(I)=U(I)+SN(1)*DU(I)/3.0
ENDDO
CALL DIF(NN,U,DU,N,M,MK)
RETURN
END SUBROUTINE RGKT2
!Runge-Kutta法右端函數值子程序
SUBROUTINE DIF(NN,U,DU,N,M,MK)
DOUBLE PRECISION U(NN),DU(NN),UUW(N*M,1)
DU(1)=1.0
DO I=2,N*M+1
DU(I)=U(I+N*M)
ENDDO
DO J=1,N*M
UUW(J,1)=U(J+1) !便于用內部函數MATMUL計算矩陣相乘
ENDDO
IF(U(1)>=0.AND.U(1)<=TD)THEN
PT=ETA1*U(1)**2+ETA2*U(1)
ELSE
PT=0
ENDIF
DO J=1,N*M
LLL1(J,1)=UUW(J,1)-AS(J,1)
ENDDO
CALL INTEGRAL2(G1,P10,F10,UUW,P10,F10,UUW,GM0,CI0,N,M,HX,HY)
CALL INTEGRAL2(G2,P0,F20,UUW,P20,F0,UUW,GM0,CI0,N,M,HX,HY)
CALL INTEGRAL2(G6,P10,F10,AS,P10,F10,AS,GM0,CI0,N,M,HX,HY)
CALL INTEGRAL2(G7,P0,F20,AS,P20,F0,AS,GM0,CI0,N,M,HX,HY)
DO J=1,N*M
L1(J,1)=E*(G1(J,1)-G2(J,1)-G6(J,1)+G7(J,1))
ENDDO
FT=MATMUL(M2N,L1)
CALL INTEGRAL2(G3,GM20,CI0,FT,P0,F20,UUW,P0,F0,N,M,HX,HY)
CALL INTEGRAL2(G4,GM0,CI20,FT,P20,F0,UUW,P0,F0,N,M,HX,HY)
CALL INTEGRAL2(G5,GM10,CI10,FT,P10,F10,UUW,P0,F0,N,M,HX,HY)
CALL INTEGRAL4(GGG1,GMB2,CI0,FT,PB,F20,UUW,PB,F0,N,M,MK,HX)
CALL INTEGRAL4(GGG2,GMB,CI20,FT,PB,F20,UUW,PB,F0,N,M,MK,HX)
CALL INTEGRAL4(GGG3,PB,F20,LLL1,PB,F20,UUW,PB,F0,N,M,MK,HX)
L2=MATMUL(M3,LLL1)
L3=MATMUL(SS5,FT)
L4=MATMUL(SS6,FT)
L5=MATMUL(SS3,LLL1)
L6=MATMUL(SS2,LLL1)
L7=MATMUL(S5,UUW)
L8=MATMUL(SS4,UUW)
DO J=1,N*M
L9(J,1)=-(L2(J,1)&
-HP*(G3(J,1)+G4(J,1)-2.0*G5(J,1))&
+(PP0+PT)*L7(J,1)&
-(EK*AK*EEK/E*(L3(J,1)-VK*L4(J,1))&
-EK*IXK*L5(J,1)&
-GK*JK*L6(J,1)&
+EK*AK/E*(GGG1(J,1)-VK*GGG2(J,1))&
-EK*AK/E/HP*(PP0+PT)*L8(J,1)&
-EK*AK*EEK*GGG3(J,1))&
-Q0*QQ0(J,1))
ENDDO
L10=MATMUL(M1N,L9)
DO I=N*M+2,2*N*M+1
DU(I)=L10(I-N*M-1,1)
ENDDO
RETURN
END SUBROUTINE DIF
END
!五次樣條函數φ5(x)
FUNCTION F5(X)
INTEGER X
DOUBLE PRECISION F5
SELECT CASE(X)
CASE(-2)
F5=1.0/120
CASE(-1)
F5=26.0/120
CASE(0)
F5=66.0/120
CASE(1)
F5=26.0/120
CASE(2)
F5=1.0/120
CASE DEFAULT
F5=0.0
END SELECT
END
!五次樣條函數φ′5(x)
FUNCTION F51(X)
INTEGER X
DOUBLE PRECISION F51
SELECT CASE(X)
CASE(-2)
F51=1.0/24
CASE(-1)
F51=5.0/12
CASE(0)
F51=0.0
CASE(1)
F51=-5.0/12
CASE(2)
F51=-1.0/24
CASE DEFAULT
F51=0.0
END SELECT
END
!五次樣條函數二階導數φ5″(x)
FUNCTION F52(X)
INTEGER X
DOUBLE PRECISION F52
SELECT CASE(X)
CASE(-2)
F52=1.0/6
CASE(-1)
F52=1.0/3
CASE(0)
F52=-1.0
CASE(1)
F52=1.0/3
CASE(2)
F52=1.0/6
CASE DEFAULT
F52=0.0
END SELECT
END
!五次樣條函數四階導數φ5⑷(x)
FUNCTION F54(X)
INTEGER X
DOUBLE PRECISION F54
SELECT CASE(X)
CASE(-2)
F54=1.0
CASE(-1)
F54=-4.0
CASE(0)
F54=6.0
CASE(1)
F54=-4.0
CASE(2)
F54=1.0
CASE DEFAULT
F54=0.0
END SELECT
END
!矩陣求逆子程序
SUBROUTINE INV(A,N,L,IS,JS)
DOUBLE PRECISION A(N,N)
INTEGER IS(N),JS(N)
L=1
DO K2=1,N
D=0.0
DO I=K2,N
DO J=K2,N
IF(ABS(A(I,J))>D)THEN
D=ABS(A(I,J))
IS(K2)=I
JS(K2)=J
ENDIF
ENDDO
ENDDO
IF(D+1.0==1.0)THEN
L=0
RETURN
ENDIF
DO J=1,N
T=A(K2,J)
A(K2,J)=A(IS(K2),J)
A(IS(K2),J)=T
ENDDO
DO I=1,N
T=A(I,K2)
A(I,K2)=A(I,JS(K2))
A(I,JS(K2))=T
ENDDO
A(K2,K2)=1/A(K2,K2)
DO J=1,N
IF(J/=K2)THEN
A(K2,J)=A(K2,J)*A(K2,K2)
ENDIF
ENDDO
DO I=1,N
IF(I/=K2)THEN
DO J=1,N
IF(J/=K2)THEN
A(I,J)=A(I,J)-A(I,K2)*A(K2,J)
ENDIF
ENDDO
ENDIF
ENDDO
DO I=1,N
IF(I/=K2)THEN
A(I,K2)=-A(I,K2)*A(K2,K2)
ENDIF
ENDDO
ENDDO
DO K2=N,1,-1
DO J=1,N
T=A(K2,J)
A(K2,J)=A(JS(K2),J)
A(JS(K2),J)=T
ENDDO
DO I=1,N
T=A(I,K2)
A(I,K2)=A(I,IS(K2))
A(I,IS(K2))=T
ENDDO
ENDDO
RETURN
END
!求板部分一重積分子程序;
SUBROUTINE INTEGRAL1(S1,FF1,FF2,N,H)
DOUBLE PRECISION S1(N,N),FF1(0:N+1,N),FF2(0:N+1,N),SS1(N,N,0:N+1),SSS1(N,N),H
INTEGER N
DO K=0,N+1
DO J=1,N
DO I=1,N
SS1(I,J,K)=FF1(K,I)*FF2(K,J)
ENDDO
ENDDO
ENDDO
SSS1=0.0
DO K=1,N
DO J=1,N
DO I=1,N
SSS1(I,J)=SSS1(I,J)+2.0*SS1(I,J,K)
ENDDO
ENDDO
ENDDO
DO I=1,N
DO J=1,N
S1(I,J)=H/2.0*(SS1(I,J,0)+SSS1(I,J)+SS1(I,J,N+1))
ENDDO
ENDDO
RETURN
END SUBROUTINE INTEGRAL1
!板部分兩個一重積分后再張量積子程序;
SUBROUTINE KRONECKER1(S2,FF3,FF4,N,HY,FF5,FF6,M,HX)
DOUBLE PRECISION FF3(0:N+1,N),FF4(0:N+1,N),FF5(0:M+1,M),FF6(0:M+1,M),FF7(N,N),FF8(M,M),&
S2(N*M,N*M),HX,HY
INTEGER N,M
CALL INTEGRAL1(FF7,FF3,FF4,N,HY)
CALL INTEGRAL1(FF8,FF5,FF6,M,HX)
DO I1=1,N
DO J1=1,N
DO I2=(I1-1)*M+1,I1*M
DO J2=(J1-1)*M+1,J1*M
S2(I2,J2)=FF7(I1,J1)*FF8(I2-(I1-1)*M,J2-(J1-1)*M)
ENDDO
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE KRONECKER1
!筋部分兩個一重積分后再張量積子程序;
SUBROUTINE KRONECKER2(SS1,FFF3,FFF4,N,MK,FFF5,FFF6,M,HX)
DOUBLE PRECISION FFF3(MK,N),FFF4(MK,N),FFF5(0:M+1,M),FFF6(0:M+1,M),FFF1(N,N),FFF2(M,M),&
SS1(N*M,N*M),HX,FFF7(N,MK)
INTEGER N,M,MK
FFF7=TRANSPOSE(FFF3)
FFF1=MATMUL(FFF7,FFF4)
CALL INTEGRAL1(FFF2,FFF5,FFF6,M,HX)
DO I1=1,N
DO J1=1,N
DO I2=(I1-1)*M+1,I1*M
DO J2=(J1-1)*M+1,J1*M
SS1(I2,J2)=FFF1(I1,J1)*FFF2(I2-(I1-1)*M,J2-(J1-1)*M)
ENDDO
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE KRONECKER2
!二重積分中的張量積子程序;
SUBROUTINE KRONECKER3(GG1,FFFF1,N,FFFF2,M)
DOUBLE PRECISION FFFF1(0:N+1,N),FFFF2(0:M+1,M),GG1((N+2)*(M+2),N*M)
INTEGER N,M
DO I1=0,N+1
DO J1=1,N
DO I2=I1*(M+2)+1,(I1+1)*(M+2)
DO J2=(J1-1)*M+1,J1*M
GG1(I2,J2)=FFFF1(I1,J1)*FFFF2(I2-I1*(M+2)-1,J2-(J1-1)*M)
ENDDO
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE KRONECKER3
!兩個張量積相乘的二重積分子程序;
SUBROUTINE INTEGRAL2(G1,FF1,FF2,A,FF3,FF4,C,FF5,FF6,N,M,HX,HY)
DOUBLE PRECISION FF1(0:N+1,N),FF2(0:M+1,M),FF3(0:N+1,N),FF4(0:M+1,M),G1(N*M,1),&
A(N*M,1),C(N*M,1),GG1((N+2)*(M+2),N*M),GG2((N+2)*(M+2),N*M),&
AA((N+2)*(M+2),1),CC((N+2)*(M+2),1),FF5(0:N+1,N),FF6(0:M+1,M),&
GG3((N+2)*(M+2),N*M),GG4((N+2)*(M+2),N*M),HX,HY
INTEGER N,M
CALL KRONECKER3(GG1,FF1,N,FF2,M)
CALL KRONECKER3(GG2,FF3,N,FF4,M)
CALL KRONECKER3(GG3,FF5,N,FF6,M)
AA=MATMUL(GG1,A)
CC=MATMUL(GG2,C)
DO J=1,N*M
DO I1=1,N
DO I2=I1*(M+2)+1,(I1+1)*(M+2)
GG3(I2,J)=2.0*GG3(I2,J)
ENDDO
ENDDO
ENDDO
DO J=1,N*M
DO I1=1,N+2
DO I2=(I1-1)*(M+2)+2,I1*(M+2)-1
GG3(I2,J)=2.0*GG3(I2,J)
ENDDO
ENDDO
ENDDO
DO J=1,N*M
DO I=1,(N+2)*(M+2)
GG4(I,J)=(HX/2.0)*(HY/2.0)*AA(I,1)*CC(I,1)*GG3(I,J)
ENDDO
ENDDO
G1=0.0
DO J=1,N*M
DO I=1,(N+2)*(M+2)
G1(J,1)=G1(J,1)+GG4(I,J)
ENDDO
ENDDO
RETURN
END SUBROUTINE INTEGRAL2
!筋部分二重積分中的張量積子程序;
SUBROUTINE KRONECKER4(G1,FF1,N,MK,FF2,M)
DOUBLE PRECISION FF1(MK,N),FF2(0:M+1,M),G1(MK*(M+2),N*M)
INTEGER N,M,MK
DO I1=1,MK
DO J1=1,N
DO I2=(I1-1)*(M+2)+1,I1*(M+2)
DO J2=(J1-1)*M+1,J1*M
G1(I2,J2)=FF1(I1,J1)*FF2(I2-(I1-1)*(M+2)-1,J2-(J1-1)*M)
ENDDO
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE KRONECKER4
!筋部分兩個張量積相乘的二重積分子程序;
SUBROUTINE INTEGRAL4(G1,FF1,FF2,A,FF3,FF4,C,FF5,FF6,N,M,MK,HX)
DOUBLE PRECISION FF1(MK,N),FF2(0:M+1,M),FF3(MK,N),FF4(0:M+1,M),G1(N*M,1),&
A(N*M,1),C(N*M,1),GG1(MK*(M+2),N*M),GG2(MK*(M+2),N*M),&
AA(MK*(M+2),1),CC(MK*(M+2),1),FF5(MK,N),FF6(0:M+1,M),&
GG3(MK*(M+2),N*M),GG4(MK*(M+2),N*M),HX
INTEGER N,M,MK
CALL KRONECKER4(GG1,FF1,N,MK,FF2,M)
CALL KRONECKER4(GG2,FF3,N,MK,FF4,M)
CALL KRONECKER4(GG3,FF5,N,MK,FF6,M)
AA=MATMUL(GG1,A)
CC=MATMUL(GG2,C)
DO J=1,N*M
DO I1=1,MK
DO I2=(I1-1)*(M+2)+2,I1*(M+2)-1
GG3(I2,J)=2.0*GG3(I2,J)
ENDDO
ENDDO
ENDDO
DO J=1,N*M
DO I=1,MK*(M+2)
GG4(I,J)=(HX/2.0)*AA(I,1)*CC(I,1)*GG3(I,J)
ENDDO
ENDDO
G1=0.0
DO J=1,N*M
DO I=1,MK*(M+2)
G1(J,1)=G1(J,1)+GG4(I,J)
ENDDO
ENDDO
RETURN
END SUBROUTINE INTEGRAL4
!載荷部分;
!求積分;
SUBROUTINE INTEGRAL6(S6,SS10,N,H)
DOUBLE PRECISION S6(N,1),SS10(0:N+1,N),SS6(N,1),H
INTEGER N
SS6=0.0
DO I=1,N
DO J=1,N
SS6(I,1)=SS6(I,1)+SS10(J,I)
ENDDO
ENDDO
DO I=1,N
S6(I,1)=H/2.0*(SS10(0,I)+2.0*SS6(I,1)+SS10(N+1,I))
ENDDO
RETURN
END SUBROUTINE INTEGRAL6
!兩個一重積分后再張量積子程序;
SUBROUTINE KRONECKER5(S2,FF3,N,HY,FF5,M,HX)
DOUBLE PRECISION FF3(0:N+1,N),FF5(0:M+1,M),FF7(N,1),FF8(M,1),&
S2(N*M,1),HX,HY
INTEGER N,M
CALL INTEGRAL6(FF7,FF3,N,HY)
CALL INTEGRAL6(FF8,FF5,M,HX)
DO I1=1,N
DO I2=(I1-1)*M+1,I1*M
S2(I2,1)=FF7(I1,1)*FF8(I2-(I1-1)*M,1)
ENDDO
ENDDO
RETURN
END SUBROUTINE KRONECKER5
!求確切點張量積子程序;
SUBROUTINE KRONECKER6(S2,FF7,N,FF8,M)
DOUBLE PRECISION FF7(N,N),FF8(M,M),S2(N*M,N*M)
INTEGER N,M
DO I1=1,N
DO J1=1,N
DO I2=(I1-1)*M+1,I1*M
DO J2=(J1-1)*M+1,J1*M
S2(I2,J2)=FF7(I1,J1)*FF8(I2-(I1-1)*M,J2-(J1-1)*M)
ENDDO
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE KRONECKER6?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -