NDSolve:Coupled PDE's, initial-boundary value problem: unreasonable “insufficient number of boundary conditions” error

I tried to NDSolve the PDE system: $$\partial_t y = x\partial_z w \quad\quad \partial_t w = \partial_z y \quad \quad \partial_z x=w$$ for $$(t,z)\in[0,1]\times[-1,0]$$ with initial conditions $$x(0,z)=x_0(z)\quad\quad w(0,z)=-sin^2(z\pi)\quad\quad y(0,z)=1$$ where $$x_0(z)$$ is the solution of $$x_0'(z)=w(0,z)\quad\quad x_0(0)=0$$ and boundary conditions $$x(t,0)=0 \quad\quad w(t,0)=w(t,-1)=0$$

From a mathematical point of view the initial and boundary conditions provided above are sufficient.

However mathematica warns that "an insufficient number of boundary conditions have been specified" and as a result "Artificial boundary effects may be present in the solution".

So there must be something wrong with my code

(*Constructing initial conditions*)
wo[z_] := -Sin[z*π]^2
yo[z_] := 1
s = NDSolve[{x'[z] == wo[z], x[0] == 0}, {x}, {z, 0, -1}]
xo[z_] := First[x[z] /. s]

(*Evolution of initial conditions towards t=1*)
equations := {D[yt[t, z], t] == xt[t,z]*D[wt[t, z], z],
D[wt[t, z], t] == D[yt[t, z], z],
D[xt[t, z], z] == wt[t, z], wt[0, z] == wo[z], yt[0, z] == yo[z],
wt[t, 0] == 0, wt[t, -1] == 0, xt[t, 0] == 0, xt[0, z] == xo[z]}
system = NDSolve[equations, {wt , yt, xt}, {z, 0, -1}, {t, 0, 1}]


I get three errors (NDSolve::pdord, NDSolve::bcart, NDSolve::ndcf), the second being "insufficient number of boundary conditions," but I cannot see what went wrong. Can anyone help?

• Welcome. What about a boundary condition for yt? Also you may want to look at the DAE tutorial. Perhaps you can find an idea there, not sure though. – user21 Nov 27 '18 at 8:47
• On the boundary $w(t,0)=0$ so $\partial_t w(t,0) = 0$ which leads to $\partial_z y(t,0)=0$. Numerically this means $y(t,0)=y(t,-\Delta x)$ where $\Delta x$ is the z-grid spacing. So an additional boundary condition would overdetermine the problem. – user61386 Nov 27 '18 at 9:01
• Version 11.3.0 produces "NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations." and "NDSolve::bcart: 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. " and " NDSolve::bcart: 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." – user64494 Nov 27 '18 at 10:54
• I agree that yt should not need a boundary condition in z. However, it is possible that Mathematica is confused by this problem and, therefore, is incorrectly asking for the boundary condition. I would try decomposing the PDEs into coupled sets of ODEs in t and solving the resulting third-order system. See the introduction to reference.wolfram.com/language/tutorial/…. As an aside, you will obtain better accuracy by solving for xo using DSolveValue instead of NDSolve. Finally, it is better not to use SetDelayed when Set will do. – bbgodfrey Nov 27 '18 at 12:16
• To continue, searching for "user:1063 do-it-yourself" will provide three examples of how to decompose PDEs into ODEs using Mathematica. I am sure that there are other examples as well on this site. – bbgodfrey Nov 27 '18 at 12:31

This computation resembles that in question 175080, in that the PDEs include spatial integrals over a dependent variable. NDSolve is not capable of solving this sort of problem as PDEs. Thus, it is necessary to perform the computation by discretizing the PDE in z, as outlined in Introduction to Method of Lines. Begin by specifying the initial conditions for w and y. (x is a matrix of w values and does not require an initial condition.)

w0[z_] = -Sin[z*Pi]^2;
y0[z_] = 1;


Next define the arrays in z for the discretized w and y.

zmin = -1; tmax = 1;
n = 100; h = -zmin/n;
W[t_] = Table[w[i][t], {i, n + 1}];
Y[t_] = Table[y[i][t], {i, n}];


Note that w and y are interleaved in z to simplify spatial derivatives. x, on the other hand, is defined symbolically as a function of w.

wave = ListCorrelate[{1, 1} h/2, W[t]];
X[t_] = Table[-Total[wave[[i ;; n]]], {i, 1, n + 1}];


(x is collocated with w.) The discretized PDEs and their initial conditions are, then,

eqw = Thread[D[W[t], t] == Join[{0}, ListCorrelate[{-1, 1}/h, Y[t]], {0}]];
eqy = Thread[D[Y[t], t] == ListCorrelate[{1, 1}/2, X[t]] ListCorrelate[{-1, 1}/h,
W[t]]];
initw = Thread[W[0] == Table[w0[(i - 1) h - 1], {i, n + 1}]];
inity = Thread[Y[0] == Table[y0[(i - 1/2) h - 1], {i, n}]];


NDSolve handles this large set of coupled ODEs without difficulty:

lines = NDSolveValue[{eqw, eqy, initw, inity}, {W[t], Y[t]}, {t, 0, tmax}],


and the results can be plotted by

wtab = Table[(i - 1) h - 1, {i, 1, n + 1}];
ParametricPlot3D[Evaluate@Thread[{wtab, t, lines[[1]]}], {t, 0, tmax},
PlotRange -> All, AxesLabel -> {"z", "t", "w"}, BoxRatios -> {2, 2, 1},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]
ytab = Table[(i - 1/2) h - 1, {i, 1, n}];
ParametricPlot3D[Evaluate@Thread[{ytab, t, lines[[2]]}], {t, 0, tmax},
PlotRange -> All, AxesLabel -> {"z", "t", "y"}, BoxRatios -> {2, 2, 1},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]


• Why Y[t_] = Table[y[i][t], {i, n}]; is used instead of Y[t_] = Table[y[i][t], {i, n+1}];? – user61386 Dec 3 '18 at 8:35
• Also I need some help with the command X[t_] = Table[-Total[wave[[i ;; n]]], {i, 1, n + 1}]; – user61386 Dec 3 '18 at 10:18
• @user61386 W and X are defined at the cell edges in z, so there are n + 1 of them. Y are defined at cell centers to improve accuracy of numerical derivatives in z, so there are n of them. X is equal to the integral from z == 0 to the cell edge for X, determined according to the trapezoidal rule. – bbgodfrey Dec 3 '18 at 16:07

There is no time-derivative of xt[t, z], which is the source of the errors (this one NDSolve::pdord being the principal one, and the other NDSolve::bcart cited by the OP apparently being a consequence of it). Differentiate the xt[] equation with respect to t and we can get @bbgodfrey's solution.

(*Evolution of initial conditions towards t=1*)
equations = {
D[yt[t, z], t] == xt[t, z]*D[wt[t, z], z],
D[wt[t, z], t] == D[yt[t, z], z],
D[D[xt[t, z], z] == wt[t, z], t],    (* differentiate xt[] equation *)
wt[0, z] == wo[z], wt[t, 0] == 0, wt[t, -1] == 0,
yt[0, z] == yo[z],
xt[0, z] == xo[z], xt[t, 0] == 0};
system = NDSolve[equations,
{wt, yt, xt}, {z, 0, -1}, {t, 0, 1},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 100, "DifferenceOrder" -> 2}},
MaxStepSize -> {0.01, 1}]  (* an upper bound on the time step size helps *)


Visualization:

ParametricPlot3D[
Evaluate@
Thread[{wt["Coordinates"] /. First@system // Last, t,
wt[t, wt["Coordinates"] /. First@system // Last]} /.
First@system], {t, 0, 1}, PlotRange -> All,
AxesLabel -> {"z", "t", "w"}, BoxRatios -> {2, 2, 1},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]


ParametricPlot3D[
Evaluate@
Thread[{yt["Coordinates"] /. First@system // Last, t,
yt[t, yt["Coordinates"] /. First@system // Last]} /.
First@system], {t, 0, 1}, PlotRange -> All,
AxesLabel -> {"z", "t", "w"}, BoxRatios -> {2, 2, 1},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]


The above elicits an error:

NDSolve::eerr: Warning: scaled local spatial error estimate...

This doesn't happen with @bbgodfrey's code, because the system is solved as a system of ODEs and not PDEs: no spatial error estimate is performed. The error can be eliminated by reducing both the spatial and time step size:

system = NDSolve[equations,
{wt, yt, xt}, {z, 0, -1}, {t, 0, 1},
Method -> {"MethodOfLines", "TemporalVariable" -> t,
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> 2}}, MaxStepSize -> 0.0025]

• Very creative (+1). – bbgodfrey Nov 30 '18 at 2:38