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. – Alex Trounev Oct 20 at 11:52
  • 1
    @alex It's {feq, fic, fbc}, isn't it? – xzczd Oct 20 at 12:30
  • @xzczd After your solving the problem, it became clear what problem is being solved. – Alex Trounev Oct 20 at 14:10
up vote 11 down vote accepted

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.

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;

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}]

Mathematica graphics

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

Mathematica graphics

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.

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 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};
(* 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;

var = Outer[#[#2, #3] &, {η, v}, grid@t, grid@y];
{barray, marray} = 
  CoefficientArrays[{Outer[#@#2 &, {ae, aeic}, {1, 2}], aebc} // Flatten, var // Flatten];
sollst = LinearSolve[marray, -barray];
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. – jsxs Oct 20 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. – xzczd Oct 20 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? – jsxs Oct 20 at 15:10
  • @jsxs Yes, just execute grid@y and observe. – xzczd Oct 20 at 15:16
  • 1
    @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 – xzczd Oct 22 at 10:42

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.