边做边学拉普拉斯，但在 Python 中得到意想不到的结果

Question

我正在尝试了解有关拉普拉斯变换的更多信息，因此我尝试在代码中实现正向和反向（梅林的逆公式）变换（使用梯形规则近似）。我希望在一个接一个地进行正向和反向时得到大致相同的信息。但是，输出值似乎与输入数据无关。

代码：

# Dependencies:
from math import ceil
from cmath import *
import numpy as np


# Constants
j = complex(0, 1)
e = exp(1).real


# Default Values
sigma_default = 0  # Real component. When 0, the result is the Fourier transform


# Forward Transform - Time Domain to Laplace Domain
def Laplace(data, is_inverse, sigma=sigma_default, frequency_stamps=None, time_stamps=None):
    # Resolve empty data scenario
    data = np.asarray(data)
    if data.size <= 1:
        return data

    # Add time data if missing
    if time_stamps is None:
        if is_inverse is False:
            time_stamps = np.arange(0, data.size)
        else:
            time_stamps = np.arange(0, data.size * 2)
    else:
        time_stamps = np.asarray(time_stamps).real
        if time_stamps.size is not data.size:
            time_stamps = np.arange(0, data.size)

    # Add frequency stamps if missing
    if frequency_stamps is None:
        if is_inverse is False:
            frequency_stamps = np.asarray(np.arange(0, ceil(data.size / 2))).real * 2 * pi  # Added forgotten constant
        else:
            frequency_stamps = np.asarray(np.arange(0, ceil(data.size))).real * 2 * pi  # Added forgotten constant
    else:
        frequency_stamps = np.asarray(frequency_stamps).real
    frequency_stamps = sigma + frequency_stamps * j

    # Create the vector of powers exp(1) is raised to. Also create the delta times / frequencies
    if is_inverse is False:
        power = -Get_Powers(time_stamps, frequency_stamps)
        delta = np.diff(time_stamps)
    else:
        power = Get_Powers(frequency_stamps, time_stamps)
        delta = np.diff(frequency_stamps)
    delta = np.concatenate([[np.average(delta)], delta])  # Ensure a start value is present

    # Perform a numerical approximation of the Laplace transform
    laplace = data * np.power(e, power) * delta

    # Trapezium rule => average 1st and last wrt zero
    laplace = laplace.transpose()  # Fixed bug in trapezium rule implementation
    laplace[[0, -1]] *= 0.5
    laplace = laplace.transpose()
    laplace = np.sum(laplace, 1)  # Integrate

    # If inverse function, then normalise and ensure the result is real
    if is_inverse is True:
        laplace *= 1 / (2 * pi * j)  # Scale
        laplace = laplace.real  # Ensure time series is real only

    # Return the result
    return laplace


# Used to derive the vector of powers exp(1) is to be raised to
def Get_Powers(values1, values2):
    # For forward Laplace, 1 = time, 2 = frequency
    # For inverse Laplace, 1 = frequency, 2 = time
    power = np.ones([values1.size, values2.size])
    power = (power * values2).transpose() * values1
    return power


if __name__ == "__main__":
    # a = [0, 1, 2, 3, 4, 5]
    a = np.arange(0, 10)
    b = Laplace(a, False)
    c = Laplace(b, True)
    print(np.asarray(a))
    print(c)

预期结果：

[0 1 2 3 4 5 6 7 8 9]
[0 1 2 3 4 5 6 7 8 9]

实际结果：

[0 1 2 3 4 5 6 7 8 9]
[162. 162. 162. 162. 162. 162. 162. 162. 162. 162.]

有什么想法我出错了吗？

编辑 1：添加拉普拉斯函数：

正向变换：

逆变换：

s的定义：

omega 在我的代码中表示为frequency_stamps。当sigma = 0时系统变成傅里叶变换。

编辑 2：修复了两个错误。问题仍然存在

Answer 1

除了原始问题中的两个错误修复之外，我还通过 Cris Luengo 的建议确定了另外 3 个错误，以调查从 Fourier Transform into the Discrete Fourier Transform 的转换。所有错误修复的摘要如下：

1. 修复了我如何实现梯形规则的错误。
2. 将frequency_stamps 缩放2*pi 以反映拉普拉斯数据的基本循环性质。
3. 再次重新调整frequency_stamps 的比例，使其仅绕一圈（也就是数据在0 -> 2*pi 范围内）。
4. 修正了一个错误，即我认为频率点只需要比时间点多一半。那是错误的。两者的数量应该相等。
5. 允许传递逆变换的初始和最终时间序列点，否则数据会损坏。

更新代码：

# Dependencies:
from cmath import *
import numpy as np


# Constants
j = complex(0, 1)
e = exp(1).real


# Default Values
sigma_default = 0.0  # Real component. When 0, the result is the Fourier transform
ends_default = np.asarray([0, 0])


# Forward Transform - Time Domain to Laplace Domain
def Laplace(data, is_inverse, sigma=sigma_default, frequency_stamps=None, time_stamps=None, ends=ends_default):
    # Resolve empty data scenario
    data = np.asarray(data)
    if data.size <= 1:
        return data

    # Add time data if missing
    if time_stamps is None:
        time_stamps = np.arange(0, data.size)  # Size doesn't change between forward and inverse
    else:
        time_stamps = np.asarray(time_stamps).real
        if time_stamps.size is not data.size:
            time_stamps = np.arange(0, data.size)

    # Add frequency stamps if missing
    if frequency_stamps is None:
        frequency_stamps = np.asarray(np.arange(0.0, data.size)).real  # Size doesn't change between forward and inverse
        frequency_stamps *= 2 * pi / np.max(frequency_stamps)  # Restrict the integral range to 0 -> 2pi
    else:
        frequency_stamps = np.asarray(frequency_stamps).real
    frequency_stamps = sigma + frequency_stamps * j

    # Create the vector of powers exp(1) is raised to. Also create the delta times / frequencies
    if is_inverse is False:
        power = -Get_Powers(time_stamps, frequency_stamps)
        delta = np.diff(time_stamps)
    else:
        power = Get_Powers(frequency_stamps, time_stamps)
        delta = np.diff(frequency_stamps)
    delta = np.concatenate([[np.average(delta)], delta])  # Ensure a start value is present

    # Perform a numerical approximation of the Laplace transform
    laplace = data * np.power(e, power) * delta
    laplace = laplace.transpose()
    laplace[[0, -1]] *= 0.5  # Trapezium rule => average 1st and last wrt zero
    laplace = laplace.transpose()
    laplace = np.sum(laplace, 1)  # Integrate

    # If inverse function, then normalise and ensure the result is real
    if is_inverse is True:
        laplace *= 1 / (2 * pi * j)  # Scale
        laplace = laplace.real  # Ensure time series is real only

        # Correct for edge cases
        laplace[0] = ends[0]
        laplace[-1] = ends[-1]

    # Return the result
    return laplace


# Used to derive the vector of powers exp(1) is to be raised to
def Get_Powers(values1, values2):
    # For forward Laplace, 1 = time, 2 = frequency
    # For inverse Laplace, 1 = frequency, 2 = time
    power = np.ones([values1.size, values2.size])
    power = (power * values2).transpose() * values1
    return power


if __name__ == "__main__":
    a = np.arange(3, 13)
    b = Laplace(a, False, sigma=0.5)
    c = Laplace(b, True, sigma=0.5, ends=np.asarray([3, 12]))
    print(np.asarray(a))
    print(c)

输出

[ 3  4  5  6  7  8  9 10 11 12]
[ 3.  4.  5.  6.  7.  8.  9. 10. 11. 12.]

感谢您的帮助！