您所在位置：主题：关于laplace逆变换的问题
jinpan2046
[专家分：0] 发布于 2009-03-13 16:00:00
程序如下：
USE MSIMSL
IMPLICIT NONE
REAL*8 Afa
REAL*8 q0D
REAL*8 k
REAL*8 Dx,TTimeD
REAL*8 RELERR,SIGMA0
INTEGER NumX, NumT, NumTT, KMAX
INTEGER I, Ix, IT, NOUT
Parameter ( KMAX=20000)
Parameter ( TTimeD=1)
Parameter ( NumX=5, NumT=1000, NumTT=NumT+1 )
REAL*8 Xc(0:NumX),Time(0:NumT),DispWs(0:NumX, 0:NumT),Monp(0:NumX, 0:NumT),Monx(0:NumX, 0:NumT),Rota(0:NumX, 0:NumT)
REAL*8 T(NumTT),FINV(NumTT)
DOUBLE COMPLEX DispL, MonpL, MonxL, RotaL
EXTERNAL DispL, MonpL, MonxL, RotaL
COMMON /Para1/Afa
COMMON /Para2/q0D
COMMON /Para3/k
COMMON /Coord/Dx
! Get output unit number
CALL UMACH (2, NOUT)
OPEN(8, FILE='Test1.txt', STATUS='UNKNOWN')
Afa =1.0
q0D =1.0
k =1.0
SIGMA0 =0.0E0
RELERR =1.0E-5
WRITE(NOUT, 3000)
WRITE(NOUT, 3020) Afa
WRITE(NOUT, 3040) k
WRITE(NOUT, 4000)
WRITE(NOUT, 4030) q0D
WRITE(NOUT, 5110)
WRITE(8, 3000)
WRITE(8, 3020) Afa
WRITE(8, 3040) k
WRITE(8, 4000)
WRITE(8, 4030) q0D
WRITE(8, 5110)
DO I=0, NumX
Xc(I)=0.0E-6+(1.0*I)/NumX
END DO
DO I=0, NumT
[color=FF0000] Time(I)=0.001E0+(I*TTimeD)/NumT[/color] END DO
DO IX=0, NumX
Dx=Xc(IX)
DO I=1, NumTT
T(I) = Time(I-1)
END DO
! Evaluate the Inverse Laplace Transformation
! CALL SINLP (DispL, NumTT, T, SIGMA0, RELERR, ERRVEC, FINV)
! CALL INLAP (DispL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
CALL DINLAP (DispL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
DispWs(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
CALL DINLAP (MonpL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
Monp(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
CALL DINLAP (MonxL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
Monx(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
CALL DINLAP (RotaL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
Rota(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
END DO
WRITE(8, 5030)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (DispWs(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (DispWs(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
WRITE(8, 5130)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (Monp(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (Monp(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
WRITE(8, 5230)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (Monx(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (Monx(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
WRITE(8, 5330)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (Rota(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (Rota(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
CLOSE(8)
3000 FORMAT(10X, 'The physical quantities '/)
3020 FORMAT(1X, 'The Frist Constant Afa = ', E10.4)
3040 FORMAT(1X, 'The second Constant K = ', E10.4)
4000 FORMAT(1X, /, 10X,'The dimensionless physical quantities '/)
4030 FORMAT(1X, 'The Load q0D = ' E10.4 )
5030 FORMAT (13X, 8x, 'Displacement of the Solid Skeleton ')
5050 FORMAT (13X, (10x, 'Coordinate of X',133X),/)
5080 FORMAT (3X, 'Time', 6x, 3x, 21(3X,E15.5),/)
5100 FORMAT (1X,E12.6, 3x, 31(3X,E15.5) )
5110 FORMAT (/,10X,'--- For Step Load -----',//)
5130 FORMAT (13X, 8x, 'Monp of the Liquid ')
5230 FORMAT (13X, 8x, 'Monx of the Solid ')
5330 FORMAT (13X, 8x, 'Rota of the Solid ')
END
DOUBLE COMPLEX FUNCTION TLoadL(S)
IMPLICIT NONE
REAL*8 q0D
DOUBLE COMPLEX S
COMMON /Para2/q0D
TLoadL=q0D/S
RETURN
END
DOUBLE COMPLEX FUNCTION C1(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C1=TLoadL(s)
RETURN
END
DOUBLE COMPLEX FUNCTION C2(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C2=-TLoadL(s)
RETURN
END
DOUBLE COMPLEX FUNCTION C3(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C3=TLoadL(s)*exp(-(Afa*s)**(1.0/2.0))/(exp((Afa*s)**(1.0/2.0))-exp(-(Afa*s)**(1.0/2.0)))
RETURN
END
DOUBLE COMPLEX FUNCTION C4(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C4=-TLoadL(s)*exp((Afa*s)**(1.0/2.0))/(exp((Afa*s)**(1.0/2.0))-exp(-(Afa*s)**(1.0/2.0)))
RETURN
END
DOUBLE COMPLEX FUNCTION C5(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C5=TLoadL(s)*(exp((Afa*s)**(1.0/2.0))*Afa*s+exp(-(Afa*s)**(1.0/2.0))*Afa*s-(Afa*s)**(1.0/2.0)*exp((Afa*s)**(1.0/2.0))+exp(-(Afa*s)**(1.0/2.0))*(Afa*s)**(1.0/2.0))/Afa/s/(exp((Afa*s)**(1.0/2.0))-exp(-(Afa*s)**(1.0/2.0)))/(Afa*s)**(1.0/2.0)
RETURN
END
DOUBLE COMPLEX FUNCTION C6(S)
IMPLICIT NONE
REAL*8 Dx,Afa,k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C6=-1.0/Afa/s*TLoadL(s)
RETURN
END
DOUBLE COMPLEX FUNCTION DispL(S)
IMPLICIT NONE
REAL*8 Dx,Afa,k
DOUBLE COMPLEX S,C1,C2,C3,C4,c5,c6
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
DispL=-(1.0/k+1.0/(Afa*s))*c2(s)*Dx-c3(s)*exp(sqrt(Afa*s)*Dx)/(Afa*s)-c4(s)*exp(-sqrt(Afa*s)*Dx)/(Afa*s)+c5(s)*Dx+c6(s)
RETURN
END
DOUBLE COMPLEX FUNCTION MonpL(S)
IMPLICIT NONE
DOUBLE COMPLEX S,C1,C2,C3,C4
REAL*8 Dx,Afa,k
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
MonpL=c1(S)+Dx*c2(S)+exp(sqrt(Afa*s)*Dx)*c3(S)+exp(-sqrt(Afa*s)*Dx)*c4(S)
RETURN
END
DOUBLE COMPLEX FUNCTION MonxL(S)
IMPLICIT NONE
DOUBLE COMPLEX S,C1,C2,C3,C4,C5,c6
REAL*8 Dx,Afa,k
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
MonxL=exp(sqrt(Afa*s)*Dx)*c3(S)+exp(-sqrt(Afa*s)*Dx)*c4(S)
RETURN
END
DOUBLE COMPLEX FUNCTION RotaL(S)
IMPLICIT NONE
DOUBLE COMPLEX S,C2,C3,C4,C5
REAL*8 Dx,Afa,k
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
RotaL=-1.0*c2(s)/(Afa*s)-exp(sqrt(Afa*s)*Dx)*c3(S)/sqrt(Afa*s)+exp(-sqrt(Afa*s)*Dx)*c4(S)/sqrt(Afa*s)+c5(s)
RETURN
END
这个程序是关于laplace逆变换的。用到了DINLAP函数。
现在的问题是，时间项（如红字所示）只能从0.1开始。如果是从0.01，0.001开始，程序就会报错。即：
The algorithm was not able to achieve the accuracy requested
within KMAX function evaluations for some T(I).
我的精度设置已经是双精度了。不知道该怎样做才行。
我不明白的是，为什么时间只能从0.1开始，不能更精确一些。因为我得出的数据是要画图的，尤其是0-1s的变化非常重要。如果只能从0.1开始，那么数据前面的变化就看不到了。
希望哪位大侠能帮我解决这个难题。
thank you！（另，时间项从0.1开始时，程序运行的非常顺畅，没有任何错误。）
USE MSIMSL
IMPLICIT NONE
REAL*8 Afa
REAL*8 q0D
REAL*8 k
REAL*8 Dx,TTimeD
REAL*8 RELERR,SIGMA0
INTEGER NumX, NumT, NumTT, KMAX
INTEGER I, Ix, IT, NOUT
Parameter ( KMAX=20000)
Parameter ( TTimeD=1)
Parameter ( NumX=5, NumT=1000, NumTT=NumT+1 )
REAL*8 Xc(0:NumX),Time(0:NumT),DispWs(0:NumX, 0:NumT),Monp(0:NumX, 0:NumT),Monx(0:NumX, 0:NumT),Rota(0:NumX, 0:NumT)
REAL*8 T(NumTT),FINV(NumTT)
DOUBLE COMPLEX DispL, MonpL, MonxL, RotaL
EXTERNAL DispL, MonpL, MonxL, RotaL
COMMON /Para1/Afa
COMMON /Para2/q0D
COMMON /Para3/k
COMMON /Coord/Dx
! Get output unit number
CALL UMACH (2, NOUT)
OPEN(8, FILE='Test1.txt', STATUS='UNKNOWN')
Afa =1.0
q0D =1.0
k =1.0
SIGMA0 =0.0E0
RELERR =1.0E-5
WRITE(NOUT, 3000)
WRITE(NOUT, 3020) Afa
WRITE(NOUT, 3040) k
WRITE(NOUT, 4000)
WRITE(NOUT, 4030) q0D
WRITE(NOUT, 5110)
WRITE(8, 3000)
WRITE(8, 3020) Afa
WRITE(8, 3040) k
WRITE(8, 4000)
WRITE(8, 4030) q0D
WRITE(8, 5110)
DO I=0, NumX
Xc(I)=0.0E-6+(1.0*I)/NumX
END DO
DO I=0, NumT
[color=FF0000] Time(I)=0.001E0+(I*TTimeD)/NumT[/color] END DO
DO IX=0, NumX
Dx=Xc(IX)
DO I=1, NumTT
T(I) = Time(I-1)
END DO
! Evaluate the Inverse Laplace Transformation
! CALL SINLP (DispL, NumTT, T, SIGMA0, RELERR, ERRVEC, FINV)
! CALL INLAP (DispL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
CALL DINLAP (DispL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
DispWs(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
CALL DINLAP (MonpL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
Monp(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
CALL DINLAP (MonxL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
Monx(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
CALL DINLAP (RotaL, NumTT, T, SIGMA0, RELERR, KMAX, FINV)
DO I=1, NumTT
Rota(IX, I-1) = FINV(I)
END DO
DO I=1, NumTT
T(I) = Time(I-1)
END DO
END DO
WRITE(8, 5030)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (DispWs(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (DispWs(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
WRITE(8, 5130)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (Monp(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (Monp(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
WRITE(8, 5230)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (Monx(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (Monx(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
WRITE(8, 5330)
WRITE(8, 5050)
WRITE(NOUT, 5080) (Xc(IX), IX=0, NumX)
WRITE(8, 5080) (Xc(IX), IX=0, NumX)
DO IT=0, NumT
WRITE(NOUT, 5100) Time(IT), (Rota(IX, IT), IX=0, NumX)
WRITE(8, 5100) Time(IT), (Rota(IX, IT), IX=0, NumX)
END DO
WRITE(8, '(//)')
CLOSE(8)
3000 FORMAT(10X, 'The physical quantities '/)
3020 FORMAT(1X, 'The Frist Constant Afa = ', E10.4)
3040 FORMAT(1X, 'The second Constant K = ', E10.4)
4000 FORMAT(1X, /, 10X,'The dimensionless physical quantities '/)
4030 FORMAT(1X, 'The Load q0D = ' E10.4 )
5030 FORMAT (13X, 8x, 'Displacement of the Solid Skeleton ')
5050 FORMAT (13X, (10x, 'Coordinate of X',133X),/)
5080 FORMAT (3X, 'Time', 6x, 3x, 21(3X,E15.5),/)
5100 FORMAT (1X,E12.6, 3x, 31(3X,E15.5) )
5110 FORMAT (/,10X,'--- For Step Load -----',//)
5130 FORMAT (13X, 8x, 'Monp of the Liquid ')
5230 FORMAT (13X, 8x, 'Monx of the Solid ')
5330 FORMAT (13X, 8x, 'Rota of the Solid ')
END
DOUBLE COMPLEX FUNCTION TLoadL(S)
IMPLICIT NONE
REAL*8 q0D
DOUBLE COMPLEX S
COMMON /Para2/q0D
TLoadL=q0D/S
RETURN
END
DOUBLE COMPLEX FUNCTION C1(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C1=TLoadL(s)
RETURN
END
DOUBLE COMPLEX FUNCTION C2(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C2=-TLoadL(s)
RETURN
END
DOUBLE COMPLEX FUNCTION C3(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C3=TLoadL(s)*exp(-(Afa*s)**(1.0/2.0))/(exp((Afa*s)**(1.0/2.0))-exp(-(Afa*s)**(1.0/2.0)))
RETURN
END
DOUBLE COMPLEX FUNCTION C4(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C4=-TLoadL(s)*exp((Afa*s)**(1.0/2.0))/(exp((Afa*s)**(1.0/2.0))-exp(-(Afa*s)**(1.0/2.0)))
RETURN
END
DOUBLE COMPLEX FUNCTION C5(S)
IMPLICIT NONE
REAL*8 Afa
REAL*8 k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C5=TLoadL(s)*(exp((Afa*s)**(1.0/2.0))*Afa*s+exp(-(Afa*s)**(1.0/2.0))*Afa*s-(Afa*s)**(1.0/2.0)*exp((Afa*s)**(1.0/2.0))+exp(-(Afa*s)**(1.0/2.0))*(Afa*s)**(1.0/2.0))/Afa/s/(exp((Afa*s)**(1.0/2.0))-exp(-(Afa*s)**(1.0/2.0)))/(Afa*s)**(1.0/2.0)
RETURN
END
DOUBLE COMPLEX FUNCTION C6(S)
IMPLICIT NONE
REAL*8 Dx,Afa,k
DOUBLE COMPLEX S,TLoadL
COMMON /Para1/Afa
COMMON /Para3/k
C6=-1.0/Afa/s*TLoadL(s)
RETURN
END
DOUBLE COMPLEX FUNCTION DispL(S)
IMPLICIT NONE
REAL*8 Dx,Afa,k
DOUBLE COMPLEX S,C1,C2,C3,C4,c5,c6
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
DispL=-(1.0/k+1.0/(Afa*s))*c2(s)*Dx-c3(s)*exp(sqrt(Afa*s)*Dx)/(Afa*s)-c4(s)*exp(-sqrt(Afa*s)*Dx)/(Afa*s)+c5(s)*Dx+c6(s)
RETURN
END
DOUBLE COMPLEX FUNCTION MonpL(S)
IMPLICIT NONE
DOUBLE COMPLEX S,C1,C2,C3,C4
REAL*8 Dx,Afa,k
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
MonpL=c1(S)+Dx*c2(S)+exp(sqrt(Afa*s)*Dx)*c3(S)+exp(-sqrt(Afa*s)*Dx)*c4(S)
RETURN
END
DOUBLE COMPLEX FUNCTION MonxL(S)
IMPLICIT NONE
DOUBLE COMPLEX S,C1,C2,C3,C4,C5,c6
REAL*8 Dx,Afa,k
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
MonxL=exp(sqrt(Afa*s)*Dx)*c3(S)+exp(-sqrt(Afa*s)*Dx)*c4(S)
RETURN
END
DOUBLE COMPLEX FUNCTION RotaL(S)
IMPLICIT NONE
DOUBLE COMPLEX S,C2,C3,C4,C5
REAL*8 Dx,Afa,k
COMMON /Para1/Afa
COMMON /Para3/k
COMMON /Coord/Dx
RotaL=-1.0*c2(s)/(Afa*s)-exp(sqrt(Afa*s)*Dx)*c3(S)/sqrt(Afa*s)+exp(-sqrt(Afa*s)*Dx)*c4(S)/sqrt(Afa*s)+c5(s)
RETURN
END
这个程序是关于laplace逆变换的。用到了DINLAP函数。
现在的问题是，时间项（如红字所示）只能从0.1开始。如果是从0.01，0.001开始，程序就会报错。即：
The algorithm was not able to achieve the accuracy requested
within KMAX function evaluations for some T(I).
我的精度设置已经是双精度了。不知道该怎样做才行。
我不明白的是，为什么时间只能从0.1开始，不能更精确一些。因为我得出的数据是要画图的，尤其是0-1s的变化非常重要。如果只能从0.1开始，那么数据前面的变化就看不到了。
希望哪位大侠能帮我解决这个难题。
thank you！（另，时间项从0.1开始时，程序运行的非常顺畅，没有任何错误。）