Trouble Plotting short term Solution on Maple

Question

I’m trying to solve a problem from Modern Mathematical Methods for Scientists and Engineers using Maple, specifically applying the Fokas method (Unified Transform) to the heat equation on the half-line. The problem setup matches Example 9.1 from the book.

As you can se bellow I tried to solve it on Mathematica and on Maple but the Maple graphs of Short-time behavior are noticeably different from the book's.

I tried the following code on maple. In Maple, the short-time solution $u(x,t)$ is noticeably different from the analytical solution shown in the book and from Mathematica's results. This discrepancy is especially apparent in plots of $u(x,t)$ for small times $ t \in [0.1, 0.6] $, where the Maple solution appears to be damped or shifted incorrectly.

What could be causing Maple's numerical integration of this oscillatory complex integral to diverge from the expected solution?

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)];

p1 := plot(x -> approx_u(x, 0.1), 0 .. 2, color = red);
p2 := plot(x -> approx_u(x, 0.2), 0 .. 2, color = green);
p3 := plot(x -> approx_u(x, 0.3), 0 .. 2, color = blue);
p4 := plot(x -> approx_u(x, 0.4), 0 .. 2, color = magenta);
p5 := plot(x -> approx_u(x, 0.5), 0 .. 2, color = RGB(1.0, 0.5, 0.));
p6 := plot(x -> approx_u(x, 0.6), 0 .. 2, color = cyan);
p0 := plot(exp(-x), x = 0 .. 2, linestyle = dash, color = black);
display([p0, p1, p2, p3, p4, p5, p6], legend = ["u(x,0)", "t=0.1", "t=0.2", "t=0.3", "t=0.4", "t=0.5", "t=0.6"], legendstyle = [location = right, font = [Helvetica, 12]], view = [0 .. 2, 0 .. 1.2], labels = ["x", "u(x,t)"], labelfont = [Helvetica, 12], axes = boxed, title = "Short-Time Behavior of u(x,t)");

Then I tried to solve the same problem on 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}]

With[{times = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6}, 
  styles = ColorData[97, "ColorList"]}, 
 Plot[Evaluate[Join[Table[uRay[x, t], {t, times}], {uRay[x, 0]}]], {x,
    0, 2}, PlotStyle -> 
   Join[styles[[;; Length[times]]], {{Black, Dashed}}], 
  PlotLegends -> 
   LineLegend[Join[Table["t = " <> ToString[t], {t, times}], {"t = 0"}],
     LegendLabel -> "Short-Time Behavior"], AxesLabel -> {"x", "u(x,t)"},
   PlotLabel -> "Short-Time Behavior of u(x,t)"]]```

dharr · Accepted Answer · 2025-08-26 03:59:09Z

I got accurate results after some changes, though I didn't attempt to optimize things. The changes were:

Did not load the VectorCalculus package. Not critical, but for fast/accurate work you don't want every * redefined, unless you need the capabilities of the package.
Removed Digits := 6;. This can degrade the accuracy and so you never want this setting. (Often hardware float accuracy applies anyway, but if you just want to change the display setting this is not the way to do it.)
Your integration limit was r = 0.. 200, which should be r = 0.. infinity. This might save time, but Maple has routines for infinite ranges, and once you make an arbitrary choice, you have no idea of the real accuracy. An alternative might be to change variables to a finite range.
I removed the integration options epsilon = 0.1e-7, method = _Gquad. There is nothing necessarily wrong here; I just wanted to let Maple make the choices.
I removed the Re(...) around this integrand. This was the most serious error, since without it the integration performs correctly but with it Maple was returning incorrect results like Float(infinity). The integrand is real as you can see by simplify(evalc(...)), which is probably what I would have done, but the expression is more complicated and does slow things down.
I removed the try catch around the integration. If you have an error, better to solve it than hide it; here the solution is the removal of the Re.
The result is very slow; I altered the plot options to calculate and display only 5 points and the calculation took about 1200 s.

Certainly some speedup is possible by altering the integration options but this is a good starting point. The values agree with the exact solution Maple finds to the PDE.

Here is the code:

with(plots);
with(ColorTools);
with(LinearAlgebra);
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; integrand := r -> (V(k1(r), x, t)*dk1(r) - V(k2(r), x, t)*dk2(r)); evalf(Int(integrand(r), r = 0 .. infinity))  end proc;
approx_u := proc(x::numeric, t::numeric) evalf(exp(-x/sqrt(2))*cos(t - x/sqrt(2)) + u1(x, t)); end proc;
p0 := plot(exp(-x), x = 0 .. 2, linestyle = dash, color = black):
p1 := plot(x -> approx_u(x, 0.1), 0 .. 2, color = red, adaptive = false, style = point, sample = [seq(0..2,numelems=5)]):
display([p0, p1], view = [0 .. 2, 0 .. 1.2]);

Collectives™ on Stack Overflow

Trouble Plotting short term Solution on Maple

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related