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

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.

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

• From your code it is not clear what problem you want to solve? Write the equation and boundary conditions. Oct 20, 2018 at 11:52
• @alex It's {feq, fic, fbc}, isn't it? Oct 20, 2018 at 12:30
• @xzczd After your solving the problem, it became clear what problem is being solved. Oct 20, 2018 at 14:10

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

• very helpful. With your pdetoode weapon, is there still any possibility to call the method of "DifferenceOrder" -> "Pseudospectral"? After testing, with difforder=2, the method fails with point up to about 90.
– lxy
Oct 20, 2018 at 14:50
• @jsxs Sadly it's impossible, because "DifferenceOrder" -> "Pseudospectral" doesn't accept symbolic input. (The reason might be symbolic calculation for "Pseudospectral" is too expensive. ) Try the new added FDM solution, this should allow you to choose a denser spatial grid. Oct 20, 2018 at 14:59
• thanks a lot! It's going to take me a while to understand the method. Btw, does the number of points include the boundaries of $y\in[-1,1]$. In other word, if points=51, we have 50 intervals?
– lxy
Oct 20, 2018 at 15:10
• @jsxs Yes, just execute grid@y and observe. Oct 20, 2018 at 15:16
• @jsxs Don't use For, try NestList, check this post for more information. As to the 2nd question, try the following and make sure you've understand the reasons for every modification: Clear[v, u1, try]; vsol[t_, y_] = v[t, y] /. v -> fsol[[2]]; u1[t_, y_] = I/\[Alpha] D[v[t, y], y] /. v -> fsol[[2]]; try[t_?NumericQ] := NIntegrate[Conjugate[{u1[t, y], vsol[t, y]}].{u1[t, y], vsol[t, y]}, {y, -1, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}]; Plot[try[t], {t, 0, Tend}] // Quiet // AbsoluteTiming Oct 22, 2018 at 10:42