0
|
2981 Views
|
9 Replies
|
1 Total Likes
View groups...
Share

Initial conditions of PDE is being ignored?

Posted 1 year ago
 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?
9 Replies
Sort By:
Posted 1 year ago
 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/1871https://mathematica.stackexchange.com/a/274214/1871To 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]]; 
Posted 1 year ago
 Hello Bryan. Do that, and I appreciate if you would tell me what they say.
Posted 1 year ago
 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!
Posted 1 year ago
 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 wasParams: 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. Attachments:
Posted 1 year ago
 What is a Mathematica Online notebook?
Posted 1 year ago
 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.
Posted 1 year ago
 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. Attachments:
Posted 1 year ago
 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.
Posted 1 year ago
 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.