Maple中半直线热方程Fokas法短时解绘图异常排查

数据老陈

2026-6-12

半直线热方程Fokas方法Maple与Mathematica计算差异原因分析

基于《Modern Mathematical Methods for Scientists and Engineers》例9.1，采用Fokas方法（统一变换）求解半直线热方程时，Maple得到的短时解（$t \in [0.1, 0.6]$）与书本解析解、Mathematica计算结果存在明显阻尼或偏移偏差，核心原因如下：

积分上限截断导致的误差：Mathematica的NIntegrate直接处理无穷积分区间{r, 0, Infinity}，会自适应判断积分收敛的有效区间；而Maple手动将积分上限截断为200，短时$t$较小时，积分核$\exp(-k^2 t)$衰减速度慢，截断上限无法覆盖全部有效积分区域，丢失的尾部积分直接导致解被过度阻尼。
积分方法适配性不足：Maple指定使用_Gquad（高斯求积）方法，该方法对高速振荡的复积分适应性差。短时$t$下，被积函数振荡剧烈，高斯求积的固定采样策略容易出现采样点不足，无法准确捕捉振荡特性，进而导致积分值偏差。Mathematica的NIntegrate会自动识别被积函数的振荡特性，调用专门的振荡积分算法，更适合这类复路径积分场景。
精度配置不匹配与误差累积：Maple设置全局精度Digits := 6，同时积分参数epsilon = 0.1e-7，两者精度要求不匹配。低全局精度会导致复运算过程中舍入误差累积，尤其是提取实部的步骤中，误差被进一步放大。Mathematica明确设置AccuracyGoal -> 6, PrecisionGoal -> 6，积分器会根据目标精度动态调整采样密度，数值稳定性更强。
实部提取时机的数值影响：Maple先计算完整复值函数，再取实部后进行积分；Mathematica则在积分内部直接对被积函数取实部。虽然数学上两者等价，但数值计算中，先取实部可以减少复运算的误差传播，尤其是振荡积分中复部的相互抵消更高效，Maple的计算顺序会导致不必要的误差积累。

对比代码

Maple代码

restart;
with(plots);
with(ColorTools);
with(LinearAlgebra);
with(Student[VectorCalculus]);
Digits := 6;
V := (k, x, t) -> -1/2*I*exp(k*x*I - k^2*t)*(1/(k - I) + 1/(k + I) - k*(1/(k^2 + I) + 1/(k^2 - I)))/Pi;
k1 := r -> r*exp(1/8*I*Pi);
k2 := r -> r*exp(7/8*I*Pi);
dk1 := r -> exp(1/8*I*Pi);
dk2 := r -> exp(7/8*I*Pi);
u1 := proc(x::numeric, t::numeric) local integrand; try integrand := r -> Re(V(k1(r), x, t)*dk1(r) - V(k2(r), x, t)*dk2(r)); evalf(Int(integrand(r), r = 0 .. 200, epsilon = 0.1e-7, method = _Gquad)); catch: 0; end try; end proc;
approx_u := proc(x::numeric, t::numeric) evalf(exp(-x/sqrt(2))*cos(t - x/sqrt(2)) + u1(x, t)); end proc;
surf := plot3d((x, t) -> approx_u(x, t), 0 .. 3, 0 .. (Student[VectorCalculus]):-`*`(2, Pi), grid = [40, 40], axes = boxed, labels = ["x", "t", "u(x,t)"], title = "Fokas Method Solution of Half-Line Heat Equation", shading = zhue);
[seq(approx_u(0, t), t = 0. .. 0.6, 0.1)];

Maple短时解绘图

Mathematica代码

(*Clear previous definitions*)ClearAll["Global`*"]

(*Define the kernel function for u1 and u2*)
V[k_, x_, 
  t_] := -I/(2 Pi) Exp[
   I k x - k^2 t]*((1/(k - I) + 1/(k + I)) - 
    k*(1/(k^2 + I) + 1/(k^2 - I)))

(* =========RAY CONTOUR METHOD=========*)(*Define rays*)
k1[r_] := r Exp[I Pi/8];
k2[r_] := r Exp[I 7 Pi/8];
dk1[r_] := Exp[I Pi/8];
dk2[r_] := Exp[I 7 Pi/8];

(*u1 from ray contour*)
u1ray[x_?NumericQ, t_?NumericQ] := 
 NIntegrate[
  Re[V[k1[r], x, t]*dk1[r] - V[k2[r], x, t]*dk2[r]], {r, 0, Infinity},
   AccuracyGoal -> 6, PrecisionGoal -> 6]

(*u2 from residue formula (exact)*)
u2ray[x_, t_] := Exp[-x/Sqrt[2]]*Cos[t - x/Sqrt[2]]

(*Full solution from ray method*)
uRay[x_?NumericQ, t_?NumericQ] := u1ray[x, t] + u2ray[x, t];

(*Plot full solution using ray contour*)Plot3D[
 uRay[x, t], {x, 0.1, 3}, {t, 0.1, 2 Pi}, PlotPoints -> 40, 
 Mesh -> None, PlotLabel -> "Full Solution via Ray Contour", 
 AxesLabel -> {"x", "t", "u(x,t)"}]

Table[uRay[0, t], {t, 0, 0.6, 0.1}]

Mathematica短时解绘图

内容的提问来源于stack exchange，提问作者Athanasios Paraskevopoulos