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

Question

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?

Answer 1

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}]

Answer 2

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]