NDSolve for unsteady Taylor-Goldstein equation, can't use MethodOfLines

I am trying to solve a relatively simple system of linear PDE in two variables (z,t). The equations are basically the Navier-stokes equations linearized around a constant background velocity profile with a source term for the temperature equation.

I have after a lot of effort gotten NDSolve to actually solve the equation, but when I try to use "MethodOfLines", I get a bunch of errors. Here is my code:

With[{U = Cos[Pi z / 1.2] - 2 Cos[Pi z / 1.2],
Q = Piecewise[{{Sin[Pi z/.4], z < .4}}, 0]},
pde1 = {
(* vorticity and theta equation *)
D[\[Omega][z, t], t]  + I k U \[Omega][z, t] -
w[z, t]  D[U, z, z]  == I k \[Theta][z, t] ,
D[\[Theta][z, t], t]  + I k U \[Theta][z, t] == -w [z, t] + Q ,
(* vorticity inversion *)

Derivative[2, 0][w][z, t] - k^2 w[z, t] == I k \[Omega][z, t],
Derivative[1, 0][w][2, t] + k w[2, t] == 0, w[0, t] == 0,
(* intial condition *)
\[Omega][z, 0] == 0, \[Theta][z, 0] ==
0, w[z, 0] == 0 };

sol = NDSolve[
pde1 /. {k -> 1}, { \[Omega], \[Theta], w}, {t, 0, 200}, {z, 0, 2},
Method -> {"PDEDiscretization" ->
{"MethodOfLines",
"TemporalVariable" -> t}}] // First]


And these are the error messages I get:

NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations. >>

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. >>

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. >>

First::nofirst:  has zero length and no first element. >>


I find this confusing for several reasons.

1. I should not need to specify an initial condition for w, since it has no time derivative and is diagnosed from the vorticity variable omega using the Laplacian solve. However, if I leave that term out, I get complaints about not specifying enough initial conditions.
2. I have specified two boundary conditions for w for a second order bvp, but it still complains. There are no other z-derivatives in my PDE. What gives?

Why am I getting these errors with "MethodOfLines" but not with Method->Automatic?

Since NDSolve seems to be getting very confused, is it possible to separate the "vorticity inversion" part of my PDE into its own NDSolve?

I really appreciate your help since I am Mathematica newb, and I am finding NDSolve to be a real pain to debug.

• Derivative[1, 0] w[2, t] + k w[2, t] == 0 should be Derivative[1, 0] [w][2, t] + k w[2, t] == 0. Please check your equations for other typos. Apr 23 '16 at 21:52
• I corrected the typo, and I get the same error. In fact, if I replace that robin boundary condition with the Dirichlet condition w[2,t] == 0, I still get the same error. Apr 23 '16 at 22:15
• If I use DirichletCondition[w[z, t] == 0, True] instead of w[0,t] ==0, w[2,t] == 0 NDSolve does return something. Now how can I get my robin condition to work? Apr 23 '16 at 22:19
• Although DirichletCondition produces a solution, it does not look like a stable solution. In fact, it is unclear to me whether NDSolve is capable of solving pde1 without preprocessing, as described in the example, Combined Elliptic-Parabolic PDE in 1D, in DAE Examples, which you might use as a basis for solving this problem. By the way, the first warning message given in the quesiton is to be expected, and the second probably arises when NDSolve attempts to eliminate w. Apr 24 '16 at 2:22
• Ok. I have been working on a manually discretized solution. I guess it was a kind of naive of me to hope for a way to solve this problem with one function call. That link is pretty helpful. It seems a little more readable then the "Method of Lines" document. Apr 24 '16 at 4:10

I'm not sure about why NDSolve fails, but managed to solve your PDE system by eliminating ω and θ first and discretizing the obtained PDE to a set of ODEs and adding proper option to NDSolve.

First, by observing the PDE system

With[{k = 1, U = Cos[(Pi z)/(12/10)] - 2 Cos[(Pi z)/(12/10)],
Q = SimplifyPWToUnitStep@Piecewise[{{Sin[(Pi z)/(4/10)], z < 4/10}}, 0]},
pde = {D[ω[z, t], t] + I k U ω[z, t] - w[z, t] D[U, z, z] ==
I k θ[z, t],
D[θ[z, t], t] + I k U θ[z, t] == -w[z, t] + Q,
Derivative[2, 0][w][z, t] - k^2 w[z, t] == I k ω[z, t]};
ic = {ω[z, 0] == 0, θ[z, 0] == 0, w[z, 0] == 0};
bc = {Derivative[1, 0][w][2, t] + k w[2, t] == 0, w[0, t] == 0}; ]


it's easy to notice ω and θ can be eliminated:

funtheta = Function[{z, t}, #] &[θ[z, t] /. First@Solve[pde[[1]], θ[z, t]]]

funomega = Function[{z, t}, #] &[ω[z, t] /. First@Solve[pde[[-1]], ω[z, t]]]

{neweq, newic,
newbc} = {pde[[2]], ic, bc} //. {θ -> funtheta, ω -> funomega} //
Simplify

(* {Sin[(5*Pi*z)/2]*(-1 + UnitStep[-(2/5) + z]) + w[z, t] +
(1/36)*(-36 + 25*Pi^2)*Cos[(5*Pi*z)/6]^2*w[z, t] +
(1/36)*I*(-72 + 25*Pi^2)*Cos[(5*Pi*z)/6]*Derivative[0, 1][w][z, t] +
Derivative[0, 2][w][z, t] + Cos[(5*Pi*z)/6]^2*
Derivative[2, 0][w][z, t] + 2*I*Cos[(5*Pi*z)/6]*
Derivative[2, 1][w][z, t] == Derivative[2, 2][w][z, t],

{w[z, 0] == Derivative[2, 0][w][z, 0],
I*(-36 + 25*Pi^2)*Cos[(5*Pi*z)/6]*w[z, 0] +
36*(Derivative[0, 1][w][z, 0] + I*Cos[(5*Pi*z)/6]*
Derivative[2, 0][w][z, 0] - Derivative[2, 1][w][z, 0]) == 0,
w[z, 0] == 0},

{w[2, t] + Derivative[1, 0][w][2, t] == 0,
w[0, t] == 0}} *)


The newic still looks a bit complicated, let's analyze it further:

First@Solve@newic
(* {Derivative[0, 1][w][z, 0] -> Derivative[2, 1][w][z, 0],
Derivative[2, 0][w][z, 0] -> 0, w[z, 0] -> 0} *)


It's not hard to notice the newic is actually equivalent to:

newic = {Derivative[0, 1][w][z, 0] == 0, w[z, 0] == 0}


Now we have a single PDE with simple i.c. and b.c., maybe NDSolve can solve it? Sadly it's not true:

lb = 0; rb = 2; tend = 200;
NDSolve[{neweq, newic, newbc}, {w}, {z, lb, rb}, {t, 0, tend}]


NDSolve::icfail

But luckily NDSolve can handle it after the PDE is discretized into a set of ODEs (the pdetoode is a general purpose function for discretizing PDE to ODE, its definition can be found here.):

xdifforder = 4; torder = 2;

points = 100;
grid = Array[# &, points, {lb, rb}];
ptoo = pdetoode[w[z, t], t, grid, xdifforder];
odevar = w /@ grid;
odeic = newic // ptoo;

odebc = With[{sf = 100}, Map[D[#, {t, torder}] + sf # &, newbc // ptoo, {2}]];
odeq = Delete[#, {{1}, {-1}}] &@(neweq // ptoo);
wsollst = NDSolveValue[{odeq, odeic, odebc}, odevar, {t, 0, tend},
Method -> {"EquationSimplification" -> "Solve"}]; // AbsoluteTiming

wsol = rebuild[wsollst, grid, 2]

ωsol = funomega /. w -> wsol

θsol = funtheta //. {w -> wsol, ω -> funomega};

Animate[Plot[wsol[x, t] // {Re@#, Im@#} & // Evaluate, {t, 0, tend},
PlotRange -> 60], {x, lb, rb}]


Remark:

1. points should be large enough.

2. SimplifyPWToUnitStep@ in the definition of Q is necessary, or NDSolve will spit out ndnum warning and fails. It's probably a bug of NDSolve.

3. The "strange" definition for odebc is necessary, if one simply forms a DAE system with the discretized b.c. i.e. use something like odebc = newbc // ptoo;, NDSolve will fail again. If you want to know more about the definition of odebc, check this obscure tutorial. (Particularly the part about Boundary Conditions. )

4. It seems to be OK to set sf to any non-negative value, I just use 100 out of habit, but do notice it's necessary to make sf > 0 when i.c. and b.c. are inconsistent.

5. It's necessary to eliminate ω and θ first, a direct discretization still fails.