Initial conditions of PDE is being ignored?

bfoo · November 10, 2022, 5:25am

I’m trying to use Mathematica to solve some PDEs numerically, and to get used to the syntax I was trying to solve some pretty trivial ones, like this:

m = 1;
k = 1;
\[Omega] = Sqrt[m^2 + k^2];
xmax = 1000;
sol[t_, x_] = \[Phi][t, x] /. NDSolve[
     {-D[\[Phi][t, x], {t, 2}] + D[\[Phi][t, x], {x, 2}] - 
        m^2*\[Phi][t, x] == 0,
      \[Phi][t, -xmax] == Exp[-I*k*xmax - I*\[Omega]*t],
      \[Phi][t, xmax] == Exp[I*k*xmax - I*\[Omega]*t],
      \[Phi][0, x] == Exp[I*k*x],
      Derivative[1, 0][\[Phi]][0, x] == -I*\[Omega]*Exp[I*k*x]
      }, \[Phi][t, x], {t, 0, 2}, {x, -xmax, xmax}][[1]];
Plot[{Re[sol[0, x]], Im[sol[0, x]]}, {x, -500, 500}]

Here xmax is just some value I set the boundaries of the solution at; the actual solution shouldn’t depend on it. I set one initial condition as \[Phi][0, x] == Exp[I*k*x] so when I plot the solution at t=0 I expect to see a plane wave with wavenumber 1, but instead I get something which is obviously different, a plane wave with wavenumber not 1:

Weirdly enough, if I lower xmax enough, it seems to plot the correct solution. Why is Mathematica not respecting the initial conditions I gave? Am I missing something crucial in the conditions I provided?

hdolhaine · November 10, 2022, 10:21pm

Hello, I think your code is quite ok and you don’t need to modify xmax etc. Look at this

m = 1;
k = 1;
\[Omega] = Sqrt[m^2 + k^2];
xmax = 1000;
sol = \[Phi] /. 
   NDSolve[{-D[\[Phi][t, x], {t, 2}] + D[\[Phi][t, x], {x, 2}] - 
        m^2*\[Phi][t, x] == 0, \[Phi][t, -xmax] == 
       Exp[-I*k*xmax - I*\[Omega]*t], \[Phi][t, xmax] == 
       Exp[I*k*xmax - I*\[Omega]*t], \[Phi][0, x] == Exp[I*k*x], 
      Derivative[1, 0][\[Phi]][0, x] == -I*\[Omega]*
        Exp[I*k*x]}, \[Phi], {t, 0, 2}, {x, -xmax, xmax}][[1, 1]];
Plot[{Re[sol[0, x]], Im[sol[0, x]]}, {x, -5, 5}]
Plot[{Re[sol[0, x]], Im[sol[0, x]]}, {x, -50, 50}]
Plot3D[Re@sol[t, x], {t, 0, 2}, {x, -5, 5}]
Plot3D[Im@sol[t, x], {t, 0, 2}, {x, -5, 5}]

But if you plot over the whole x-range you will get a bad pic.

bfoo · November 11, 2022, 1:05am

Could you possibly attach a few images of what you get? I don’t seem to be getting any plots with the correct wavenumber, regardless of what x range I plot over.

hdolhaine · November 11, 2022, 5:07pm

Hello Bryan, I attach a notebook produced on my system (Version 7). It contains some images, perhaps you would like to look at them before trying to run the notebook.

Untitled-1.nb (1.29 MB)

bfoo · November 11, 2022, 6:49pm

Thank you. I wonder if this is a bug with my version of Mathematica, as your images look correct but if I try to run the notebook, I get different images.

EDIT: Could you try running your code on a Mathematica Online notebook? I get the same issue there as in my desktop version.

hdolhaine · November 12, 2022, 2:49pm

What is a Mathematica Online notebook?

hdolhaine · November 16, 2022, 4:03pm

Hello Bryan,

I forwarded my notebook to a friend who has a more recent version of Mathematica, and he too, doesn’t get the same pictures as I do. Seems there indeed is something wrong.
However, he pointed out that your differential equation in the NDSolve-statement might be overdetermined. In principle, you have twice the statement containing xmax…

You should delete one of them and try your code again.

I have no idea why my version (Mathematica 7) works properly and also with one of these conditions deleted.

Your original form was

Params:

m = 1;
k = 1;
\[Omega] = Sqrt[m^2 + k^2];
xmax = 1000;

giving the function f0

sol = \[Phi] /. NDSolve[
     {-D[\[Phi][t, x], {t, 2}] + D[\[Phi][t, x], {x, 2}] - 
        m^2*\[Phi][t, x] == 0, \[Phi][t, -xmax] == 
       Exp[-I*k*xmax - I*\[Omega]*t],
      \[Phi][t, xmax] == Exp[I*k*xmax - I*\[Omega]*t],
      \[Phi][0, x] == Exp[I*k*x],
      Derivative[1, 0][\[Phi]][0, x] == -I*\[Omega]*Exp[I*k*x]},
     \[Phi], {t, 0, 2}, {x, -xmax, xmax}][[1, 1]];
f0 = sol

In shorter form (one of the xmax-conditions deleted)

Clear[sol];
sol = \[Phi] /. NDSolve[
     {-D[\[Phi][t, x], {t, 2}] + D[\[Phi][t, x], {x, 2}] - 
        m^2*\[Phi][t, x] == 0,
      \[Phi][t, xmax] == Exp[I*k*xmax - I*\[Omega]*t],
      \[Phi][0, x] == Exp[I*k*x],
      Derivative[1, 0][\[Phi]][0, x] == -I*\[Omega]*Exp[I*k*x]},
     \[Phi], {t, 0, 2}, {x, -xmax, xmax}][[1, 1]];
f1 = sol

with function f1 as a result.

However, there exists a symbolic solution.
Your differential equation reads

deq[f_] := -D[f, t, t] + D[f, x, x] - m^2 f

I found a symbolic solution ff

ff = (1/2 (E^(-I Sqrt[a^2 + m^2] t) + E^(I Sqrt[a^2 + m^2] t)) C[1] + 
     1/2 I (E^(-I Sqrt[a^2 + m^2] t) - E^(I Sqrt[a^2 + m^2] t)) C[
       2]) (E^(I a x) K[1] + E^(-I a x) K[2]);

Indeed

deq[ff] // FullSimplify
(* 0 *)

Adjusting the constants, you get function ff1, which fulfills (almost) all your conditions:

ff /. t -> 0
rr1 = {K[1] -> 1, K[2] -> 0, C[1] -> 1, 
  C[2] -> -I \[Omega] /Sqrt[k^2 + m^2], a -> k}
ff1 = ff /. rr1
ff1 /. t -> 0
\[Phi][0, x] == Exp[I*k*x]

and

D[ff1, t] /. t -> 0
-I*\[Omega]*Exp[I*k*x]

Finally, fix k and m

ff1 /. {k -> 1, m -> 1}

and you will see that the two numerical solutions obtained above and the symbolic solution are identical. Plot all three functions together for different tmes and regions (try the Im-parts as well)

Manipulate[
 Plot[
  {Re@f0[tt, xx], Re@f1[tt, xx], Re@ff1 /. {t -> tt, x -> xx}},
  {xx, -10, 10}, 
  PlotStyle -> {{Dashing[.1], Thickness[.01], Black}, 
    Blue, {Dashing[.05], Thick, Red}}
  ],
 {tt, 0, 2}
 ]

I take this as a sign that NDSolve (on my system) works properly. What is happening with yours?

I attach the latest notebook as well if you consider sending it to Wolfram support.

Untitled-wcom.nb (26.5 KB)

bfoo · November 24, 2022, 7:39pm

Hello, sorry for the late response. I ran your notebook; the numerical and analytic solutions do not match up in my version of Mathematica. I sent a note to Wolfram Support about this as I’m pretty sure it’s a Mathematica bug. Thanks for your help!

hdolhaine · November 25, 2022, 9:24pm

Hello Bryan. Do that, and I appreciate if you would tell me what they say.

xzczd · November 26, 2022, 6:12am

The i.c. is not ignored. It’s because NDSolve hasn’t used enough grid points for spatial discretization. Hans doesn’t observe this phenomenon in version 7 because there’s a silent change in prior error estimation of NDSolve in version 9. This issue is discussed in:

https://mathematica.stackexchange.com/a/262173/1871

https://mathematica.stackexchange.com/a/274214/1871

To resolve the problem in higher versions, force a denser grid:

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}};
sol[t_, x_] = ϕ[t, x] /. 
   NDSolve[{-D[ϕ[t, x], {t, 2}] + D[ϕ[t, x], {x, 2}] - m^2*ϕ[t, x] == 
       0, ϕ[t, -xmax] == Exp[-I*k*xmax - I*ω*t], ϕ[t, xmax] == 
       Exp[I*k*xmax - I*ω*t], ϕ[0, x] == Exp[I*k*x], 
      Derivative[1, 0][ϕ][0, x] == -I*ω*Exp[I*k*x]}, ϕ[t, x], {t, 0,
       2}, {x, -xmax, xmax}, Method -> mol[10^4, 4]][[1]];

Topic		Replies	Views
NDSolve Overwrites Initial Condition in Linearized PDE Questions astronomy , mathematica	0	3259	October 1, 2014
Error with Numerical Solution of PDE Questions numerical-computation	2	3583	September 7, 2015
Can Mathematica modify lightly the inital conditions of NDSolve?? Questions mathematica , calculus , equation-solving , numerical-computation	4	7446	January 15, 2015
DSolve and initial condition Questions equation-solving , symbolic-computations	5	6342	September 12, 2020
Get solution of a PDE with NDSolve? Questions equation-solving , numerical-computation	2	7172	September 10, 2019

Initial conditions of PDE is being ignored?

Related topics