SymPy逆拉普拉斯变换绘图空轴无报错问题排查与解决

问题场景 #

在Colab环境中使用SymPy执行逆拉普拉斯变换计算，调用plot函数绘制变换结果时，出现空坐标轴但无任何报错信息。相关代码如下：

import numpy as np
from sympy import *
import sympy as sp
from sympy.physics.units.quantities import Quantity
import sys
sys.set_int_max_str_digits(0)

%config InlineBackend.figure_format = 'svg'

def laplace(f):
    s = symbols('s')
    t = symbols ('t',real=True)
    return sp.laplace_transform(f, t, s, noconds=True)

def ilaplace(F):
    s = symbols('s')
    t = symbols ('t',real=True)
    return sp.inverse_laplace_transform(F, s, t,noconds=True).evalf().nsimplify()

s,n = symbols('s,n')

Rc= sp.Rational(1)
Rp= sp.Rational(1,2)
Cw= sp.Rational(1,5)
Cs= sp.Rational(1,200)
Cl= sp.Rational(1,5)
T= sp.Rational(4)
P1= sp.Rational(5)
Q= sp.Rational(1,3)

t = symbols('t',real = True) 

Ct=(1/Cl + 1/Cw + 1/Cs)**(-1)
A=1 / (Rp*Ct)
B=Rc
C=1/Cs + Rc/(Rp*Ct)
D=(1/(Rp*Cs))*(1/Cl+1/Cw)

r1=(-C+sqrt(C**2-4*B*D))/(2*B)
r2=(-C-sqrt(C**2-4*B*D))/(2*B)

P=Sum(laplace(P1*(sp.Heaviside(t-n*T)-sp.Heaviside(t-T*(n+Q)))),[n,0,1])
G=(s+A)/((s-r1)*(s-r2))

V=(G*P)
v=ilaplace(V)
v=v
v

grafV1=plot(v,(t,0,50),xlabel="t",ylabel="V(t)",legends=True,line_color="green")

问题原因 #

1. 逆变换函数的过早数值化：ilaplace函数中使用.evalf().nsimplify()将符号表达式强制转换为数值形式，破坏了SymPy绘图所需的符号结构，导致无法正确解析表达式生成数据点。
2. Sum对象未展开：P是未求值的Sum对象，逆拉普拉斯变换无法完全处理该结构，最终得到的v表达式包含未解析的求和项，绘图时无法计算有效数值。
3. Heaviside函数的分段处理：逆变换后生成的分段函数（由Heaviside函数转换而来）未被SymPy的plot函数正确识别，采样时遗漏了有效区间的数据。

解决方法 #

1. 修改逆拉普拉斯变换函数 #

移除.evalf().nsimplify()，保留纯符号表达式，避免过早数值化破坏结构：

def ilaplace(F):
    s = symbols('s')
    t = symbols ('t',real=True)
    return sp.inverse_laplace_transform(F, s, t, noconds=True)

2. 显式展开Sum求和 #

在定义P时调用.doit()展开求和，得到具体的表达式而非Sum对象：

P = Sum(laplace(P1*(sp.Heaviside(t-n*T)-sp.Heaviside(t-T*(n+Q)))), [n,0,1]).doit()

3. 整理分段函数（可选） #

对逆变换后的v进行展开或分段整理，确保表达式结构清晰：

v = ilaplace(V)
v = v.expand()  # 展开表达式
v = sp.piecewise_fold(v)  # 整理分段函数结构

4. 调整绘图采样参数 #

关闭自适应采样并增加采样点数量，确保覆盖所有关键区间：

grafV1 = plot(v, (t, 0, 50), xlabel="t", ylabel="V(t)", legends=True, line_color="green", adaptive=False, nb_of_points=1000)

完整修正代码 #

import numpy as np
from sympy import *
import sympy as sp
from sympy.physics.units.quantities import Quantity
import sys
sys.set_int_max_str_digits(0)

%config InlineBackend.figure_format = 'svg'

def laplace(f):
    s = symbols('s')
    t = symbols ('t',real=True)
    return sp.laplace_transform(f, t, s, noconds=True)

def ilaplace(F):
    s = symbols('s')
    t = symbols ('t',real=True)
    return sp.inverse_laplace_transform(F, s, t, noconds=True)

s,n = symbols('s,n')

Rc= sp.Rational(1)
Rp= sp.Rational(1,2)
Cw= sp.Rational(1,5)
Cs= sp.Rational(1,200)
Cl= sp.Rational(1,5)
T= sp.Rational(4)
P1= sp.Rational(5)
Q= sp.Rational(1,3)

t = symbols('t',real = True) 

Ct=(1/Cl + 1/Cw + 1/Cs)**(-1)
A=1 / (Rp*Ct)
B=Rc
C=1/Cs + Rc/(Rp*Ct)
D=(1/(Rp*Cs))*(1/Cl+1/Cw)

r1=(-C+sqrt(C**2-4*B*D))/(2*B)
r2=(-C-sqrt(C**2-4*B*D))/(2*B)

# 展开Sum求和
P=Sum(laplace(P1*(sp.Heaviside(t-n*T)-sp.Heaviside(t-T*(n+Q)))),[n,0,1]).doit()
G=(s+A)/((s-r1)*(s-r2))

V=(G*P)
v=ilaplace(V)
v = v.expand()
v = sp.piecewise_fold(v)

# 调整绘图参数
grafV1=plot(v,(t,0,50),xlabel="t",ylabel="V(t)",legends=True,line_color="green", adaptive=False, nb_of_points=1000)

内容的提问来源于stack exchange，提问作者Simón Pedro Durá