NDSolve of a stiff system of ODE's

Question

I am trying to solve numerically the following set of differential equations

eqs = { 
D L^2 c'[y] - L my[y] == 0, 
D L^2 my''[y] - L (Qyy'[y] + c'[y]/3) - my[y] == 0, 
D L^2 Qyy''[y] - (4/15) L my'[y] - Qyy[y] == 0,
c[-1] - L my[-1] - 1/(3 D) == 0,
D L^2 my'[-1] - L (Qyy[-1] + (1 + c[-1])/3) == 0, 
D L^2 my'[1] - L (Qyy[1] + (1 + c[1])/3) == 0,
D L^2 Qyy'[-1] - (4/15) L my[-1] == 0, 
D L^2 Qyy'[1] - (4/15) L my[1] == 0
}

What I am actually interested is the solution to above set of equations for the parameters D and L being from the interval (0,1). For not too small values of D, and L, one can just simply use NDSolve without any further problems to obtain an accurate solution.

For small values of D and L, for example D = 1/2, and L = 1/10, the solution is very stiff close to the boundaries y = - 1, and 1, so that the solution is problematic

pars = {D -> 1/2, L -> 1/10};
sol = NDSolve[eqs /. pars, {c, my, Qyy}, {y, -1, 1}];
Plot[c[y] /. sol, {y, -1, 1}, PlotRange -> All]

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly.
NDSolve::berr: The scaled boundary value residual error of 2.9891304910534367`*^10 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

First approach

My first idea was to solve this problem using StiffnessSwitching method, but this does not automatically solve the problem

pars = {D -> 1/2, L -> 1/10};
sol = NDSolve[eqs /. pars, {c, my, Qyy}, {y, -1, 1}, 
   Method -> "StiffnessSwitching"];
Plot[c[y] /. sol, {y, -1, 1}, PlotRange -> All]

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly.
NDSolve::berr: The scaled boundary value residual error of 1.631995982197461`*^8 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

I assume the problem is related to the mesh that StiffnessSwitching method is using. I could in principle set MaxStepSize -> 0.001, but I guess this makes the whole StiffnessSwitching method useless if the MaxStepSize is smaller than the step size used in StiffnessSwitching. Is there a different way of approaching this problem? Is there a way of telling mathematica to use a finer mesh, but to preserve the uneven distribution of points convenient for dealing with stiff functions?

Second approach

I also tried to solve this problem using spectral methods. If I understand the documentation correctly (https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) the following code should solve the equations on the Chebyshev grid, but the result does not seem to differ from the original approach that was done without specifying any Method at all

pars = {D -> 1/2, L -> 1/10};
points = 80;
sol = NDSolve[eqs /. pars, {c, my, Qyy}, {y, -1, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "DifferenceOrder" -> "Pseudospectral", "MaxPoints" -> points, 
       "MinPoints" -> points}}];
Plot[c[y] /. sol, {y, -1, 1}, PlotRange -> All]

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly.
NDSolve::berr: The scaled boundary value residual error of 2.9891304910534367`*^10 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

Is there something I am missing in here? Even varying the number of points does not seem to have any effect on the result.

General coding advice: (Avoid starting your symbols with capitals)[https://mathematica.stackexchange.com/a/18395/4999], especially single capital letters, and more especially single capital letters like D that are protected system commands. — Michael E2, Sep 02 '21 at 14:37
Try WorkingPrecision -> 32. Raising working precision is the rich person's solution to bad conditioning. You might also add PrecisionGoal -> 8 to keep the precision goal the same. It generally will be set to half the working precision. Raising the precision goal tends to make the conditioning seem worse, but the practical reason is that you probably don't need a more precise result. — Michael E2, Sep 02 '21 at 14:47
Thanks, changing WorkingPrecision, and PrecisionGoal indeed helps to obtain a sensible result for these set of parameters. However, if I want to solve the problem for D = 1/10, and L = 1/10, the same problem appears. — user35780, Sep 02 '21 at 15:31
WorkingPrecision -> 100, PrecisionGoal -> 8 worked for me in that case. It does suggest it may keep getting more difficult as D approaches zero. — Michael E2, Sep 03 '21 at 01:11
Yes this seem to work, thanks! Would there also be a way of doing this using spectral methods (as described in my second approach)? — user35780, Sep 03 '21 at 14:12
I don't know of a built-in way in Mathematica for ODEs, even though they have it for the spatial component of a PDE. The Chebfun package in MATLAB can do it (see section 7.8). — Michael E2, Sep 03 '21 at 16:45

bbgodfrey · Answer 1 · 2021-09-05T00:07:34.137

Because this system of equations is linear with constant coefficients, it can be solves symbolically.

sol = DSolveValue[eqs, {c[y], qyy[y], my[y]}, {y, -1, 1}];

With

LeafCount@sol
(* 73666 *)

sol is too lengthy to reproduce here. It can, however, be plotted without difficulty. For instance,

Plot[Evaluate[sol /. {d -> .1, l -> .1}], {y, -1, 1}, PlotRange -> {-1, 1}, 
  ImageSize -> Large,  AxesLabel -> {y, "c,qyy,my"}, LabelStyle -> {15, Bold, Black}]

Values near the endpoints are exponentially large, literally.

sol /. {d -> .1, l -> .1, y -> -1}
(* {-3.60777*10^11, -1.18112*10^11, 4.46677*10^11} *)
sol /. {d -> .1, l -> .1, y -> 1}
(* {-4.12317*10^11, -1.18112*10^11, -4.98216*10^11} *)

It is understandable that NDSolve has progressively more difficulty solving these equations as d and l approach zero.

NDSolve of a stiff system of ODE's

1 Answers1