如何用Plots.jl基于solve()输出数据绘制绝对值平方曲面图

最小工作示例（MWE） #

using DifferentialEquations

N₁=5
γ=-1 # 非线性项强度参数

g₁(u,p,t) = 1*im*(p[1]*Tridiagonal(fill(1,N₁), fill(-2,N₁+1), fill(1,N₁))*u + γ*u.^3)
vals=0.3*ones(N₁+1)
u0  = map(Complex, vals)

tspan = (0.0,1.0)
prob = ODEProblem(g₁,u0,tspan, [0.6, 1])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

需求 #

基于上述ODE求解得到的sol数据，绘制绝对值平方的三维曲面图：横轴对应离散空间网格点，纵轴对应时间，竖轴为每个时空点上解的绝对值平方，效果为典型的时空演化三维曲面。

尝试过程与问题 #

尝试使用Plots.jl的surface()函数绘图，执行代码：

using Plots
surface(range(-5,5), sol.t, abs.(sol))

但触发错误，核心原因：

• sol是ODESolution类型，直接对其调用abs.()无法生成曲面图所需的二维数据结构
• 空间范围range(-5,5)默认生成100个点，与解的空间维度（此处为6个点）不匹配，维度无法对应

解决方法 #

需先将sol转换为符合曲面图要求的二维矩阵，计算绝对值平方，同时保证空间轴的长度匹配：

步骤1：转换数据结构 #

将sol.u中每个时间点的复数向量堆叠为矩阵，每行对应一个时间点，每列对应一个空间网格点：

# 将sol.u的向量堆叠为空间×时间的矩阵，再转置为时间×空间的结构
sol_matrix = hcat(sol.u...)'
# 计算每个元素的绝对值平方（abs2比abs再平方更高效）
abs2_matrix = abs2.(sol_matrix)

步骤2：定义匹配的空间轴 #

确保空间范围的点数与解的空间维度一致（此处为N₁+1=6个点）：

x = range(-5, 5, length=N₁+1)

步骤3：绘制曲面图 #

surface(x, sol.t, abs2_matrix, 
        xlabel="空间", ylabel="时间", zlabel="|u|²",
        title="绝对值平方的时空演化曲面")

完整可运行代码 #

using DifferentialEquations, Plots

N₁=5
γ=-1 # 非线性项强度参数

g₁(u,p,t) = 1*im*(p[1]*Tridiagonal(fill(1,N₁), fill(-2,N₁+1), fill(1,N₁))*u + γ*u.^3)
vals=0.3*ones(N₁+1)
u0  = map(Complex, vals)

tspan = (0.0,1.0)
prob = ODEProblem(g₁,u0,tspan, [0.6, 1])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

# 数据处理
sol_matrix = hcat(sol.u...)'
abs2_matrix = abs2.(sol_matrix)
x = range(-5, 5, length=N₁+1)

# 绘图
surface(x, sol.t, abs2_matrix, 
        xlabel="空间", ylabel="时间", zlabel="|u|²",
        title="绝对值平方的时空演化曲面")

内容的提问来源于stack exchange，提问作者KZ-Spectra