SymPy Piecewise积分生成3D线框图与手动计算结果不符排查

问题场景 #

在求解一维热方程的傅里叶级数近似时，使用SymPy的Piecewise函数定义分段被积函数并积分，生成的3D线框图结果与手动分段计算积分的结果不符，手动分段版本为正确结果。代码中1到10的循环用于替代傅里叶级数的无穷求和（级数会收敛）。

错误代码（Piecewise版本） #

import numpy as np
import sympy as sm
from sympy import *
from spb import *

x = sm.symbols("x")
t = sm.symbols("t")
n = sm.symbols("n", integer=True)

L = 20
D = 0.475
f = (S(2)/L)*sin(n*np.pi*x/20)*Piecewise((x, (0 <= x) & (x <= 10)), (20-x, (10 < x) & (x <= 20)))
print('The function u(x,0) : ')
print('')
sm.pretty_print(f)
print('')
print('')
print(piecewise_fold(f))

fint = integrate(f, (x, 0, 20))
g = fint*exp(-(n**2)*(np.pi**2)*D*t/400).nsimplify()
print('')
print('')
sm.pretty_print(fint)
print('')
print('')

s3 = 0
for c in range(10):
    s3 += g.subs({n:c})

print('')
print('The function u(x,t) : ')
print('')
sm.pretty_print(s3)

plot3d(
    s3, (x, 0, 20), (t, 0, 10), {"alpha": 0}, # hide the surface
    wireframe=True, wf_n1=20, wf_n2=10,
    wf_rendering_kw={"color": "tab:blue"}, # optional step to customize the wireframe lines
    backend=MB, zlabel="$u(x,t)$", title="One Dimensional Heat Equation"
)

正确代码（手动分段版本） #

import numpy as np
import sympy as sm
from sympy import *
from spb import *

x = sm.symbols("x")
t = sm.symbols("t")
n = sm.symbols("n", integer=True)

L = 20
f1 = (2/L)*x*sin(n*np.pi*x/20)
f2 = (2/L)*(20-x)*sin(n*np.pi*x/20)
fint1 = sm.integrate(f1,(x,0,10))
fint2 = sm.integrate(f2,(x,10,20))

D = 0.475
g = (fint1+fint2)*sin(n*np.pi*x/20)*exp(-(n**2)*(np.pi**2)*D*t/400).nsimplify()
s = 0
for c in range(10):
    s += g.subs({n:c})
    print(s)

print('')
print('The function u(x,t) : ')
print('')
sm.pretty_print(s)
print('')
print('')

plot3d(
    s, (x, 0, 20), (t, 0, 10), {"alpha": 0}, # hide the surface
    wireframe=True, wf_n1=20, wf_n2=10,
    wf_rendering_kw={"color": "tab:blue"}, # optional step to customize the wireframe lines
    backend=MB, zlabel="$u(x,t)$", title="One Dimensional Heat Equation"
)

问题根源 #

1. 傅里叶级数项结构错误：
   错误代码中，将sin(nπx/20)包含在了被积分的函数f中，积分得到的fint已经是傅里叶系数，但后续构建级数项时，没有再乘以sin(nπx/20)，导致最终的级数缺少了关键的空间正弦项，与热方程的解结构不符。

2. 手动分段的正确性：
   手动分段代码中，先计算分段函数与sin(nπx/20)乘积的积分（得到傅里叶系数fint1+fint2），再将系数乘以sin(nπx/20)和时间衰减项，这完全符合一维热方程傅里叶正弦级数解的形式：
   $$u(x,t) = \sum_{n=1}^{\infty} b_n \sin\left(\frac{n\pi x}{L}\right) e^{-\frac{n2\pi^2 D t}{L^2}}$$
   其中$b_n = \frac{2}{L} \int_0^L u(x,0) \sin\left(\frac{n\pi x}{L}\right) dx$，而$u(x,0)$是分段函数。

修正后的Piecewise版本代码 #

import numpy as np
import sympy as sm
from sympy import *
from spb import *

x = sm.symbols("x")
t = sm.symbols("t")
n = sm.symbols("n", integer=True)

L = 20
D = 0.475
# 定义初始分段函数u(x,0)
u0 = Piecewise((x, (0 <= x) & (x <= 10)), (20-x, (10 < x) & (x <= 20)))
# 计算傅里叶系数：b_n = (2/L)∫₀^L u0 * sin(nπx/L) dx
b_n = (S(2)/L) * integrate(u0 * sin(n*np.pi*x/L), (x, 0, L))
# 构建级数项：b_n * sin(nπx/L) * 时间衰减项
g = b_n * sin(n*np.pi*x/L) * exp(-(n**2)*(np.pi**2)*D*t/(L**2)).nsimplify()

s3 = 0
for c in range(1, 10):  # 注意从1开始，n=0时项为0，可跳过
    s3 += g.subs({n:c})

print('The function u(x,t) : ')
sm.pretty_print(s3)

plot3d(
    s3, (x, 0, 20), (t, 0, 10), {"alpha": 0},
    wireframe=True, wf_n1=20, wf_n2=10,
    wf_rendering_kw={"color": "tab:blue"},
    backend=MB, zlabel="$u(x,t)$", title="One Dimensional Heat Equation"
)

内容的提问来源于stack exchange，提问作者Freya the Goddess