Return to Revisions

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.

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

So, let's turn to the method shown here i.e. discretize the system ourselves and make NDSolve use the ODE solver to solve the problem. 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.

Clear[η0, v0, α, β]
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}, Map[sf # + D[#, t] &, fbc, {3}]];
var = Outer[#[#2] &, {η, v}, grid];

sollst = NDSolveValue[{ode /@ {1, 2}, odeic, odebc}, var, {t, 0, Tend}];

(* A more advanced but faster approach: *)
(*
lhs = D[Flatten[var][t] // Through, t];
{barray, marray} = CoefficientArrays[{ode /@ {1, 2}, odebc} // Flatten, lhs];
rhs = LinearSolve[marray, -barray]; // AbsoluteTiming

sollst = NDSolveValue[{lhs == rhs, odeic}, var, {t, 0, Tend}];
 *)

rebuild = ListInterpolation[
    Developer`ToPackedArray@#["ValuesOnGrid"] & /@ # // 
     Transpose, {#[[1]]["Coordinates"][[1]], grid}] &;

fsol = rebuild /@ sollst;

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, but notice this method will probably fail if their values are too large, because it'll be harder to transform the system to the standard form as they become larger.

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.