# Mathematica and IBVP mixing temporal and spatial derivatives

The following PDE's-system-solving code

 cn = 10^-1; zmin = -1*cn; tmax = 1*cn;
IBVP = NDSolve[{D[w[t, z], t] ==
D[q[t, z], z] + w[t, z]*D[w[t, z], z],
D[q[t, z], t] == D[w[t, z], z] + w[t, z]*D[q[t, z], z],
w[0, z] == 0, w[t, 0] == 0, q[0, z] == 1,
q[t, 0] == Cos[t]^2}, {w, q}, {t, 0, tmax}, {z, 0, zmin},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 100}}];
Plot3D[w[t, z] /. IBVP, {t, 0, tmax}, {z, 0, zmin}]


is just fine for Mathematica. However when the PDE $$\partial_zB(t,z)=0$$ is added to the system along with trivial initial-boundary conditions $$B(0,z)=0\qquad B(t,0)=0$$ i.e. the code becomes

cn = 10^-1; zmin = -1*cn; tmax = 1*cn;
IBVP2 = NDSolve[{D[w[t, z], t] ==
D[q[t, z], z] + w[t, z]*D[w[t, z], z],
D[q[t, z], t] == D[w[t, z], z] + w[t, z]*D[q[t, z], z],
D[B[t, z], z] == 0, w[0, z] == 0, w[t, 0] == 0, q[0, z] == 1,
q[t, 0] == Cos[t]^2, B[0, z] == 0, B[t, 0] == 0}, {w, q, B}, {t,
0, tmax}, {z, 0, zmin},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 100}}];
Plot3D[w[t, z] /. IBVP2, {t, 0, tmax}, {z, 0, zmin}]


Mathematica displays the following error

Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable z. Artificial boundary effects may be present in the solution. >>

It is conceivable that one can solve the IBVP for $$B(t,z)$$ seperately and get $$B(t,z)=0$$.

However it may be the case that D[B[t, z], z] has a dependence on w[t, z] or q[t, z] with D[w[t, z], t] having also a dependance on B[t, z] in which case the PDE's are coupled and can not be integrated indpendently from one another.

So it is technically important to find out what is wrong with IBVP2.

I would appreciate any help on the above.

PS: My real concern is to tackle problems where the PDE's are coupled and mix temporal and spatial derivatives for example $$\partial_tw=e^B\cdot\partial_zq+w\cdot\partial_zw\qquad\partial_tq=e^B\cdot\partial_zw+w\cdot\partial_zq\qquad\partial_zB=w\cdot e^q$$ so I need an answer that can be generalised in this case.

I am not sure whether Μathematica can tackle problems of this kind in general. However it is hard to believe that it crushes for this -apparently- trivial case.

• In short, NDSolve can't handle PDE that doesn't involve pure derivative of time well (at least now) . Similar problems: mathematica.stackexchange.com/q/184281/1871 mathematica.stackexchange.com/q/163923/1871 mathematica.stackexchange.com/q/183745/1871 mathematica.stackexchange.com/q/118194/1871 mathematica.stackexchange.com/q/133731/1871 Jan 9 '19 at 8:53
• However there are some exceptions when NDSolve can handle PDE that doesn't involve pure derivative of time. Have you seen this post? Jan 9 '19 at 9:10
• Please notice Alex has modified i.c. and b.c. drastically in that example, which is usually impossible in real cases. Jan 9 '19 at 9:11
• How can such a trivial example confuse Mathematica? There must be some way to fix the code. Jan 9 '19 at 12:16
• In your PS, I think your current BCs and ICs will always make everything zero except for q at the boundary (both B and w are zero initially making it impossible to evolve). Can you set up a case with your desired BCs and ICs? Jan 9 '19 at 15:16

There is a simple modification that fits the occasion. You just have to remove the condition B[t,0] == 0 . Then the message appears, but the two solutions for w[t,z], q[t,z] although not identical, but not very different.

cn = 10^-1; zmin = -1*cn; tmax = 1*cn;
IBVP = NDSolve[{D[w[t, z], t] ==
D[q[t, z], z] + w[t, z]*D[w[t, z], z],
D[q[t, z], t] == D[w[t, z], z] + w[t, z]*D[q[t, z], z],
w[0, z] == 0, w[t, 0] == 0, q[0, z] == 1,
q[t, 0] == Cos[t]^2}, {w, q}, {t, 0, tmax}, {z, 0, zmin},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 100}}];

IBVP2 = NDSolve[{D[w[t, z], t] ==
D[q[t, z], z] + w[t, z]*D[w[t, z], z],
D[q[t, z], t] == D[w[t, z], z] + w[t, z]*D[q[t, z], z],
D[B[t, z], z] == 0, w[0, z] == 0, w[t, 0] == 0, q[0, z] == 1,
q[t, 0] == Cos[t]^2, B[0, z] == 0}, {w, q, B}, {t, 0, tmax}, {z,
0, zmin},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 100}}];

{Plot3D[w[t, z] /. IBVP, {t, 0, tmax}, {z, 0, zmin}, PlotRange -> All,
Mesh -> None, ColorFunction -> Hue, AxesLabel -> {"t", "z", "w"},
PlotLabel -> "IBVP"],
Plot3D[w[t, z] /. IBVP2, {t, 0, tmax}, {z, 0, zmin},
PlotRange -> All, Mesh -> None, ColorFunction -> Hue,
AxesLabel -> {"t", "z", "w"}, PlotLabel -> "IBVP2"],
Plot[{w[t, -t] /. IBVP, w[t, -t] /. IBVP2}, {t, 0, tmax},
PlotLegends -> {"IBVP", "IBVP2"}]} Consider a solution with b.c. B[t, 0] == 0

IBVP1 = NDSolve[{D[w[t, z], t] ==
D[q[t, z], z] + w[t, z]*D[w[t, z], z],
D[q[t, z], t] == D[w[t, z], z] + w[t, z]*D[q[t, z], z],
D[B[t, z], z] == 0, w[0, z] == 0, w[t, 0] == 0, q[0, z] == 1,
q[t, 0] == Cos[t]^2, B[0, z] == 0, B[t, 0] == 0}, {w, q, B}, {t,
0, tmax}, {z, 0, zmin},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 100}}];
{Plot3D[q[t, z] /. IBVP, {t, 0, tmax}, {z, 0, zmin}, PlotRange -> All,
Mesh -> None, ColorFunction -> Hue, AxesLabel -> {"t", "z", "q"},
PlotLabel -> "IBVP"],
Plot3D[q[t, z] /. IBVP1, {t, 0, tmax}, {z, 0, zmin},
PlotRange -> All, Mesh -> None, ColorFunction -> Hue,
AxesLabel -> {"t", "z", "q"}, PlotLabel -> "IBVP1"],
Plot[{q[t, -t] /. IBVP, q[t, -t] /. IBVP1}, {t, 0, tmax},
PlotLegends -> {"IBVP", "IBVP1"}]}


In this case, the two solutions IBVP, IBVP1 for q[t,z] diverge. But the two solutions IBVP2, IBVP1 for B[t,z] are the same B[t,z]=0. I found a method for solving a system of equations in the general case. We use the explicit Euler in time and the standard solver for z. We solve the example given by the author as "My real concern":

zmin = -1/10; t0 = 1/20; tmax = 63*t0;
W[z_] := 0
Q[z_] := 1
B[z_] := 0
Do[{W[t], Q[t], B[t]} =
NDSolveValue[{(w[z] - W[t - t0][z])/t0 ==
Exp[b[z]]*D[q[z], z] + w[z]*D[w[z], z], (q[z] - Q[t - t0][z])/
t0 == Exp[b[z]]*D[w[z], z] + w[z]*D[q[z], z],
D[b[z], z] == w[z]*Exp[q[z]], w == 0, q == Cos[t]^2,
b == 0}, {w, q, b}, {z, 0, zmin}, Method -> "BDF"], {t, t0,
tmax, t0}];
{ListPlot3D[
Flatten[Table[{t, z, W[t][z]}, {t, 0, tmax, t0}, {z, 0, zmin,
zmin/50}], 1], Mesh -> None, InterpolationOrder -> 3,
ColorFunction -> "SouthwestColors", AxesLabel -> {"t", "z", "w"}],
ListPlot3D[
Flatten[Table[{t, z, Q[t][z]}, {t, 0, tmax, t0}, {z, 0, zmin,
zmin/50}], 1], Mesh -> None, InterpolationOrder -> 3,
ColorFunction -> "SouthwestColors", AxesLabel -> {"t", "z", "q"}],
ListPlot3D[
Flatten[Table[{t, z, B[t][z]}, {t, 0, tmax, t0}, {z, 0, zmin,
zmin/50}], 1], Mesh -> None, InterpolationOrder -> 3,
ColorFunction -> "SouthwestColors", AxesLabel -> {"t", "z", "B"}]} • My real concern is the case when the PDE's are coupled. In that case not giving b.c. for $B$ would be mathematically problematic. So what alternatives do I have? Jan 9 '19 at 13:09
• Your real case should be studied but not by this example. In this example, Mathematica adds artificial boundary conditions, so you need to remove b.c. B[t,0] == 0. Jan 9 '19 at 13:49
• The solution for B in IBVP2 is almost meaningless, isn't it? Without the b.c. in z direction, solution for D[B[t, z], z] == 0 won't be determined. Jan 9 '19 at 13:58
• Check that the numerical solution B = 0 is obtained with and without b.c. B[t,0] == 0. See update. Jan 9 '19 at 14:25
• Not that different. Notice the variable B is not coupled with q and w, so it should be uniquely determined by {D[B[t, z], z] == 0, B[t, 0] == (* something *)}. NDSolve happens to find B[t, z] == 0 in IBVP2, but it's no more than a coincidence. Jan 9 '19 at 15:08