01拉普拉斯变换定义

1.变换公式

2. 常见函数的Laplace变换

02 Laplace数值逆变换

根据拉普拉斯逆变换的公式，可以看到，f(t)可以变成如下的公式。

#!/usr/local/bin/python
# -*- coding: gbk -*-
#============================================================
# TEST1.PY                     -- by Dr. ZhuoQing 2020-07-10
#
# Note:
#============================================================

from headm import *

#################################################################
#   Function InvLap(t,omega,sigma,nint), numerically inverts a  #
#   Laplace transform F(s) into f(t) using the Fast Fourier     #
#   Transform (FFT) algorithm for a specific time "t", an       #
#   upper frequency limit "omega", a real parameter "sigma"     #
#   and the number of integration intervals "nint" .            #
#                                                               #
#   Function F(s) is defined in separate as Fs(s) (see code     #
#   below). Fs(s) has to be changed accordingly everytime the   #
#   user wants to invert a different function.                  #
#                                                               #
#   I suggest to use omega>100 and nint=50*omega. The higher    #
#   the values for omega, the more accurate the results will be #
#   in general, but at the expense of longer processing times.  #
#                                                               #
#   Sigma is a real number which must be a little bigger than   #
#   the real part of rightmost pole of the function F(s). For   #
#   example, F(s) = 1/s + 1/(s-2) + 1/(s+1) has poles for s=0,  #
#   s=2 and s=-1. Hence, sigma must be made equal to, say,      #
#   2.05 so as to keep all poles at the left of this value.     #
#   The analytical inverse for this simple function is          #
#   f(t) = 1 + exp(-t) + exp(2t). For t=1.25, omega=200,        #
#   nint=10000 and sigma=2.05, the numerical inversion yields   #
#   f(1.25) ~= 13.456844516, or -0.09% away from the actual     #
#   analytical result, 13.468998757 (results truncated to 9     #
#   decimal places). This is the example used in this code.     #
#                                                               #
#   Creator: Fausto Arinos de Almeida Barbuto (Calgary, Canada) #
#   Date: May 18, 2002                                          #
#   E-mail: fausto_barbuto@yahoo.ca                             #
#                                                               #
#   Reference:                                                  #
#   Huddleston, T. and Byrne, P: "Numerical Inversion of        #
#   Laplace Transforms", University of South Alabama, April     #
#   1999 (found at http://www.eng.usouthal.edu/huddleston/      #
#   SoftwareSupport/Download/Inversion99.doc)                   #
#                                                               #
#   Usage: invoke InvLap(t,omega,sigma,nint), for t>0.          #
#                                                               #
#################################################################
#   We need cmath because F(s) is a function operating on the
#   complex argument s = a + bj
from math import ceil
from cmath import *

#   *** Driver InvLap function  ***
def InvLap(t,omega,sigma,nint):
#   Sanity check on some parameters.
    omega = ceil(omega)
    nint = ceil(nint)

    if omega <= 0:
       omega = 200

    if nint <= 0:
        nint = 10000

    return (trapezoid(t,omega,sigma,nint))

#   *** Function trapezoid computes the numerical inversion. ***
def trapezoid(t,omega,sigma,nint):
    sum = 0.0
    delta = float(omega)/nint
    wi = 0.0

#   The for-loop below computes the FFT Inversion Algorithm.
#   It is in fact the trapezoidal rule for numerical integration.
    for i in range(1,(nint+1)):
        witi = complex(0,wi*t)

        wf = wi + delta
        wfti = complex(0,wf*t)

        fi = (exp(witi)*Fs(complex(sigma,wi))).real
        ff = (exp(wfti)*Fs(complex(sigma,wf))).real
        sum = sum + 0.5*(wf-wi)*(fi+ff)
        wi = wf

    return ((sum*exp(sigma*t)/pi).real)

#   *** The Laplace function F(s) is defined here.  ***
def Fs(s):
#    return (1.0/s + 1.0/(s+1.0) + 1.0/(s-2.0))
    return 1.0/s

#   Function InvLap(t,omega,sigma,nint) is invoked.
printf(InvLap(1.25, 200, 3.05, 10000))

#------------------------------------------------------------
#        END OF FILE : TEST1.PY
#============================================================

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
 

03一些常见函数数值逆变换结果

• f(t)=sin(t)
  

def Fs(s):
	return 1/(s*s+1)

  

 

  
1
  
2
 

• f(t)=t
  
  

  ^{▲ 函数t的Laplace数值逆变换结果}

• f(t)=exp(-t)

※ 结论

本文测试了 一个简易版本的数值Laplace逆变换的程序。

这个算法需要人工给出积分的上下限，积分限的实部等。如果参数给定不合适，计算出来的数值误差较大。

文章来源: zhuoqing.blog.csdn.net，作者：卓晴，版权归原作者所有，如需转载，请联系作者。

原文链接：zhuoqing.blog.csdn.net/article/details/107241738