Why does NDSolve fail to solve the PDEs and spit out mconly warning?

Question

I try to solve two coupled PDEs with NDSolve using the following code:

1. Set two operators:

   op1[y_, α_, β_] = ((α^2 + β^2)*# - D[#, {y, 2}]) &;
   op2[y_, α_, β_] = (op1[y, α, β]@ op1[y, α, β]@#) &;
   

2. Set the parameters:

   α = 1; β = 0.5; m = 300; Tend = 20; nx = 201;
   

3. Set the equation, boundary and initial conditions (BCs and ICs):

   With[{η = η[t, y], v = v[t, y], U = 1 - y^2},
    feq = {D[η, t] + I*α*U*η + 
    1/m*op1[y, α, β][η] == -I*β*D[U, y]*v,
   op1[y, α, β][D[v, t]] + 
    I*α*U*op1[y, α, β][v] + 
    I*α*D[U, {y, 2}]*v + 1/m*op2[y, α, β][v] == 0};
   fic = {η == η0[y], v == v0[y]} /. t -> 0;
   fbc = {{v == 0, η == 0, D[v, y] == 0} /. y -> -1,
   {v == 0, η == 0, D[v, y] == 0} /. y -> 1};]
   

4. Take random (smooth) functions as the ICs (In my real problem, I need random values of the functions to initiate the numerical integration in NDSolve):

   SeedRandom[1];
   η0[y_] = BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 10], {0.}], SplineClosed -> False][(y + 1)/2];
   SeedRandom[2];
   v0[y_] = BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 10], {0.}], SplineClosed -> False][(y + 1)/2];
   

5. Solve the system:

   fsol = NDSolve[{feq, fic, fbc}, {η, v}, {y, -1, 1}, {t, 0, Tend},
   Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", 
     "MinPoints" -> nx, "MaxPoints" -> nx, "DifferenceOrder" -> "Pseudospectral"}}]
   

The code has been run in Mathematica V9.0, which gives the following warnings:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

This error is understandable because the ICs are two random functions, although I managed to satisfy $v(\pm1,t)=\eta(\pm,t)=0$, the remaining $\frac{\partial v}{\partial y}(\pm1,t)=0$ have not in general. But I guess this warning is not a big deal.

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. >>

Why NDSolve is unable to deal with complex values? I do need complex numbers in the numerical solutions.

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.

This error should be the most important one but I cannot understand this issue.

6. If the PDEs could be solved, I want to plot:

   u1[t, y] = I/α*D[v[t, y], y] /. fsol[[1]];
   u3[t, y] = I/α*η[t, y] /. fsol[[1]];
   

Please help me. Any suggestion is highly appreciated!

Answer 1

Explanation for the Warning

Notice the description for the warning NDSolve::mconly is

NDSolve::mconly: For the method IDA, only machine real code is available.

What's method IDA? It's a method for solving DAE. (Please search in the document for more information. )

Why is NDSolve calling a DAE solver? Because NDSolve discretizes the PDE system based on method of lines ("MethodOfLines") and the resulting system is a DAE system or an ODE system.

Why does NDSolve choose to discretize the PDE system to a DAE system rather than an ODE system? Because your PDE system involves terms like $\frac{\partial ^2v}{\partial t\partial y}$ and current implementation of method of lines in NDSolve isn't clever enough to transform the discretized system to the standard form required by the ODE solver. (For more information about the standard form, check this post. )

Why does the DAE solver fail? Because generally the DAE solver of Mathematica is weaker than the ODE solver, at least up to now.

Solutions

Partly NDSolve-based solution

One possible method to resolve the problem is to discretize the system ourselves and help NDSolve to use the ODE solver as I've done here. I'll use pdetoode for the task.

Notice I've modified the definition of η0 and v0 a little, because pdetoode can only handle Listable function.

(* Tested in v9.0.1 *)
Clear[η0, v0, α, β]
op1[y_, α, β] = ((α^2 + β^2)*# - D[#, {y, 2}]) &;
op2[y_, α, β] = (op1[y, α, β]@ op1[y, α, β]@#) &;
SetAttributes[#, Listable] & /@ {η0, v0};
SeedRandom[1];
η0[y_?NumericQ] = 
  BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 10], {0.}], 
    SplineClosed -> False][(y + 1)/2];
SeedRandom[2];
v0[y_?NumericQ] = 
  BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 10], {0.}], 
    SplineClosed -> False][(y + 1)/2];
α = 1; β = 0.5; m = 300; Tend = 20; nx = 201;
With[{η = η[t, y], v = v[t, y], U = 1 - y^2}, 
   feq = {D[η, t] + I α U η + op1[y, α, β][η]/m == (-I) β D[U, y] v, 
          op1[y, α, β][D[v, t]] + 
            I α U op1[y, α, β][v] + I α D[U, {y, 2}] v + op2[y, α, β][v]/m == 0}; 
 fic = {η == η0[y], v == v0[y]} /. t -> 0; 
 fbc = {{v == 0, η == 0, D[v, y] == 0} /. y -> -1, 
        {v == 0, η == 0, D[v, y] == 0} /. y -> 1}; ]
domain = {-1, 1};
difforder = 2;
points = 50;
grid = Array[# &, points, domain];
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[{η, v}[t, y], t, grid, difforder];
delone = #[[2 ;; -2]] &;
deltwo = #[[3 ;; -3]] &;
ode@1 = delone@ptoofunc@feq[[1]];
ode@2 = deltwo@ptoofunc@feq[[2]];
odeic = ptoofunc@fic;
odebc = ptoofunc@With[{sf = 1}, diffbc[t, sf]@fbc];
var = Outer[#[#2] &, {η, v}, grid];
sollst = NDSolveValue[{ode /@ {1, 2}, odeic, odebc}, var, {t, 0, Tend}, 
  Method -> {"EquationSimplification" -> "MassMatrix"}];
fsol = rebuild[#, grid]& /@ sollst;
Clear[u1, u3]

u1[t_, y_] = I/α D[v[t, y], y] /. v -> fsol[[2]];
u3[t_, y_] = I/α η[t, y] /. η -> fsol[[1]];
Plot3D[u1[t, y] // Abs // Evaluate, {t, 0, Tend}, {y, -1, 1}]

Plot3D[u3[t, y] // Abs // Evaluate, {t, 0, Tend}, {y, -1, 1}]

You may try larger difforder or points.

If you want to make the solution fit the b.c. better, choose a larger sf. As to the meaning of sf, check this post.

Purely FDM-based solution

As already mentioned, transforming the system to the standard form required by ODE solver can be time consuming. (difforder = 2; points = 100 is already challenging for the method above. ) So MassMatrix method shown above turns out to be very efficient for the problem. Still, it's not a bad idea to leave NDSolve alone and turn to pure finite difference method (FDM) as I've done here. I'll use pdetoae for the task:

(* Definitions for feq etc. are the same as above. *)
Clear[domain, points, grid];
domain@t = {0, Tend}; domain@y = {-1, 1};
points@t = 100; points@y = 100;
difforder = 4;
(grid@# = Array[# &, points@#, domain@#]) & /@ {t, y};
var = Outer[#[#2, #3] &, {η, v}, grid@t, grid@y];
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[{η, v}[t, y], grid /@ {t, y}, difforder];
delone = #[[2 ;; -2]] &;
deltwo = #[[3 ;; -3]] &;
ae@1 = delone /@ Rest@ptoafunc@feq[[1]];
ae@2 = deltwo /@ Rest@ptoafunc@feq[[2]];
aeic@1 = delone@ptoafunc@fic[[1]];
aeic@2 = deltwo@ptoafunc@fic[[2]];
aebc = ptoafunc@fbc;
{barray, marray} = 
  CoefficientArrays[{Outer[#@#2 &, {ae, aeic}, {1, 2}], aebc} // Flatten, var // Flatten];
sollst = LinearSolve[marray, -barray];

Remark

If you have difficulty in understanding the usage of Rest, delone and deltwo, the following is an alternative that doesn't require you to remove equations from the system:

fullsys = ptoafunc@{feq, fic, fbc};
{barray, marray} = CoefficientArrays[fullsys // Flatten, var // Flatten];
sollst = LeastSquares[marray, -barray, Method -> Direct]; // AbsoluteTiming

Notice this approach is slower.

solfunclst = 
  ListInterpolation[#, grid /@ {t, y}] & /@ ArrayReshape[sollst, {2, points@t, points@y}];
Clear[u1, u3]
u1[t_, y_] = I/α D[v[t, y], y] /. v -> solfunclst[[2]];
u3[t_, y_] = I/α η[t, y] /. η -> solfunclst[[1]];

The resulting pictures look the same as above so I'd like to omit them here.