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

Question

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.

Answer 1

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 = Simplify`PWToUnitStep@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[(5Piz)/2](-1 + UnitStep[-(2/5) + z]) + w[z, t] + 
       (1/36)(-36 + 25Pi^2)Cos[(5Piz)/6]^2w[z, t] + 
       (1/36)I(-72 + 25Pi^2)Cos[(5Piz)/6]Derivative[0, 1][w][z, t] + 
       Derivative[0, 2][w][z, t] + Cos[(5Piz)/6]^2*
         Derivative[2, 0][w][z, t] + 2ICos[(5Piz)/6]*
         Derivative[2, 1][w][z, t] == Derivative[2, 2][w][z, t],
{w[z, 0] == Derivative[2, 0][w][z, 0], 
     I(-36 + 25Pi^2)Cos[(5Piz)/6]w[z, 0] + 
         36(Derivative[0, 1][w][z, 0] + ICos[(5Piz)/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}, diffbc[{t, torder}, sf]@newbc // ptoo];
odeq = Delete[#, {{1}, {-1}}] &@(neweq // ptoo);
wsollst = NDSolveValue[{odeq, odeic, odebc}, odevar, {t, 0, tend}, 
    Method -> {"EquationSimplification" -> "MassMatrix"}]; // AbsoluteTiming
(* Timing in v9.0.1: 11.5454 second )
( Timing in v12.3.1: 18.9494 second *)
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. Simplify`PWToUnitStep@ 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.