# Numerical solution to PDE seems to contradict initial values

I am trying to solve this PDE system:

EPS = NDSolveValue[{
D[e[z, t], z] == 0,
D[p[z, t], t] == -p[z, t] + I s[z, t],
D[s[z, t], t] == I p[z, t],
e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0
}, {e, p, s}, {z, 0, 1}, {t, 0, 1}];


Then I seek values for e and p at different z and t=0:

{EPS[[1]][0, 0], EPS[[2]][0, 0]}

{EPS[[1]][1, 0], EPS[[2]][1, 0]}


The result contradicts the condition p[z, 0] == 0:

{1, I}
{1, I}


Furthermore, if I ONLY change the order of the equations

EPS = NDSolveValue[{
D[p[z, t], t] == -p[z, t] + I s[z, t],
D[s[z, t], t] == I p[z, t],
D[e[z, t], z] == 0,
e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0
}, {e, p, s}, {z, 0, 1}, {t, 0, 1}];

{EPS[[1]][0, 0], EPS[[2]][0, 0]}

{EPS[[1]][1, 0], EPS[[2]][1, 0]}


The result changes, and is wrong again

{1, 0.999909 I}
{1.76683*10^10, - 1.09719*10^6 I}


Does anyone have an idea what is happening?

UPDATE: The original problem arised from this fully coupled equation:

NDSolveValue[{
D[e[z, t], z] == I p[z, t],
D[p[z, t], t] == -p[z, t] + I e[z, t] + I g[t] s[z, t],
D[s[z, t], t] == I Conjugate[g[t]] p[z, t],
s[z, 0] == p[z, 0] == 0, e[0, t] == f[t]
}, {e, p, s}, {z, 0, 1}, {t, 0, 1}]


Here g[t] and f[t] are simple functions, such as const, Gaussian, etc.

• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. It is recommended that you browse the Markdown help and How to copy code from Mathematica so it looks good on this site. Once read the edit window help button ? gives quick reminders of that material. To remove from the screen, click ? again. Sep 1, 2016 at 14:25
• Thanks for your note, now I hope it reads better. Sep 1, 2016 at 15:32

The first thing I tried was to check if this is actually time integrated and not treated as a pure spatial PDE. Just calling a variable t does not tell NDSolve that something is to be considered time dependent.

EPS = NDSolve[{
D[e[z, t], z] == 0,
D[p[z, t], t] == -p[z, t] + I s[z, t],
D[s[z, t], t] == I p[z, t],
e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0}, {e, p, s}, {z, 0,
1}, {t, 0, 1}, Method -> "MethodOfLines"][[1]]

NDSolve::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.


OK, so the PDE is treated as a pure spatial problem and it used the FEM for that:

EPS = NDSolve[{
D[e[z, t], z] == 0,
D[p[z, t], t] == -p[z, t] + I s[z, t],
D[s[z, t], t] == I p[z, t],
e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0}, {e, p, s}, {z, 0,
1}, {t, 0, 1}][[1]];
(e /. EPS)["ElementMesh"]

NDSolveFEMElementMesh[{{0., 1.}, {0.,
1.}}, {NDSolveFEMQuadElement["<" 400 ">"]}]


Now, we look at the equation once more and when we eliminate the constant derivative condition we get:

EPS = NDSolve[{
(*D[e[z,t],z]\[Equal]0,*)
D[p[z, t], t] == -p[z, t] + I s[z, t],
D[s[z, t], t] == I p[z, t],
e[0, t] == 1, p[z, 0] == 0, s[z, 0] == 0}, {e, p, s}, {z, 0,
1}, {t, 0, 1}][[1]]
{e -> Function[{z, t}, 1],
p -> InterpolatingFunction[{{0., 1.}, {0., 1.}}, <>],
s -> InterpolatingFunction[{{0., 1.}, {0., 1.}}, <>]}

{(e /. EPS)[0, 0], (p /. EPS)[0, 0]}
{(e /. EPS)[1, 0], (p /. EPS)[1, 0]}
{1, 0.}
{1, 0.}


This time a time integration did happen and for the spatial discretization you can use either FEM or TGP by givin the option (if you want): Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid"}} or "FiniteElement".

In the pure spatial problem a slightly different PDE is solved. Now, the reordering is an unfortunate issue that is documented and there is currently no way to reliably autodetect that and warn about it. For an explanation and workarounds see this tutorial in the section about Ordering of Dependent Variable Names.

• Thanks, I see now it works as expected. However, my original pde is more complex and the equations are coupled. The same weird result is obtained there too. How should I handle it? NDSolveValue[{ D[e[z, t], z] == I p[z, t], D[p[z, t], t] == -p[z, t] + I e[z, t] + I g[t] s[z, t], D[s[z, t], t] == I Conjugate[g[t]] p[z, t], s[z, 0] == p[z, 0] == 0, e[0, t] == f[t] }, {s, p, e}, {z, 0, 1}, {t, 0, 1}] Sep 1, 2016 at 15:36
• @SvetoslavIvanov, you could edit your question with an update and add your new PDE, also include g and f. Sep 1, 2016 at 16:12