# Stiff System of ODEs: NDSolve Convergence Issues

I have tried to use NDSolve for solving 4th order coupled ODEs problems (see the attached codes). The MMA solver fails to solve this problem.

ClearAll["Global*"]
L = 10;
ode1 = y''[t] - 0.01 y''''[t] == 0;
ic11 = y == 0;
ic12 = y'' == 0;
ic13 = y'[L] == 0;
ic14 = x[L]*x[L]*(y'[L] - 0.01*y'''[L]) == 1/20;
ode2 = -10 (4.7169 (1.12 - x[t]) + 0.2120 x''[t]) +
424000 x[t] (y'[t]^2 + 0.010 y''[t]^2) == 0;
ic21 = x' == 0;
ic22 = x'[L] == 0;
sn = NDSolveValue[{ode1, ode2, ic11, ic12, ic13, ic14, ic21,
ic22}, {x[t], y[t]}, {t, 0, L},
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {x == 1}}];


Note that the initial value of x must equal 1.12.

How can I set the "Shooting" method in MMA for solving Stiff System of ODEs? Namely, how to define the initial values for the "Shooting" method in MMA.

Update Version 01_2020.

Now we are trying to verify the method proposed by @bbgodfrey

Such test procedure looks like that:

ic4 has changed to

ic14 = x[L]x[L](y'[L] - 0.01*y'''[L]) == 0;

The input file:

L = 10;
ode1 = y''[t] - 0.01 y''''[t] == 0;
ic11 = y == 0;
ic12 = y'' == 0;
ic13 = y'[L] == 0;
ic14 = x[L]*x[L]*(y'[L] - 0.01*y'''[L]) == 0;
ode2 = -10 (4.7169 (1.12 - x[t]) + 0.2120 x''[t]) +
424000 x[t] (y'[t]^2 + 0.010 y''[t]^2) == 0;
ic21 = x' == 0;
ic22 = x'[L] == 0;
sy = (DSolve[{ode1, ic11, ic12, ic13}, y, t] // Flatten) /.
C -> c Exp[-100]
ode2x = Simplify[ode2 /. sy];

ic14x = Collect[ic14 /. sy, x, Simplify];
sn = NDSolveValue[{ode2x /. c -> c[t], ic14x /. c -> c[L], ic21, ic22,
c'[t] == 0}, {x[t], c}, t,
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {x == -1/2, c == I/8}}] //
Flatten;
sn // Last
Plot[Evaluate@ReIm@First@sn, {t, 0, L}, ImageSize -> Large,
AxesLabel -> {t, x}, LabelStyle -> {15, Bold, Black}]
Plot[Evaluate@ReIm@Last[y /. sy /. c -> Last[sn]], {t, 0, L},
ImageSize -> Large, AxesLabel -> {t, y},
LabelStyle -> {15, Bold, Black}]


the output: should be y ==0 and x = const 1.12 for ic4 (new one):

ic14 = x[L]x[L](y'[L] - 0.01*y'''[L]) == 0;

The simulated results: Obviously, x is not const, oscillation can be observed.

• I recommend that you solve for y , using any boundary conditions that depend on y only, and then insert that result into the ODE for x. Also, be sure to replace any C[_] constants in the expression for y by constants without indices before that substitution. – bbgodfrey Aug 12 at 22:21
• I added a monitor: Print[Dynamic@{foo, Clock[Infinity]}]; sn = NDSolveValue[..., StepMonitor :> (foo = t)]. It shows that shooting for the correct initial conditions takes a long time. You say the problem is stiffness — maybe you know that theoretically — but if you could give better starting initial conditions, you might solve it. [Added: @bbgodfrey's idea seems better; didn't see it.] – Michael E2 Aug 12 at 22:23
• @bbgodfrey in the second step (your method): Solve::inex: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. – ABCDEMMM Aug 12 at 23:04
• @ABCDEMMM Requiring that x == 1 as a boundary condition would overspecify this ode system. Note, however, that "StartingInitialConditions" -> {x == 1} does not guarantee that the solution for x is equal to 1 at t = 0. – bbgodfrey Aug 12 at 23:59

This system is similar in some respects to question 228267, but much more challenging. It can be solved as follows. First, as a good practice, Rationalize all equations in the question. Then, solve for y, which can be done symbolically.

sy = (DSolve[{ode1, ic11, ic12, ic13}, y, t] // Flatten) /. C -> c Exp[-100]
(* {y -> Function[{t}, 1/100 E^(-100 - 10 t)(-((E^100 c)/E^100) + (E^(100 + 20 t) c)/E^100
- (10 E^(10 t) t c)/E^100 - (10 E^(200 + 10 t) t c)/E^100)]} *)


Because c is an arbitrary constant, it can be renormalized as desired, here to assure that c is of order unity in magnitude. Applying this result to ode2 and ic14 yields

ode2x = Simplify[ode2 /. sy]
(* (47169 x[t])/1000 + 4240 E^(-20 (20 + t)) (2 E^200 + E^(20 t) + 2 E^(40 (5 + t)) -
2 E^(10 (10 + t)) + 2 E^(20 (10 + t)) - 2 E^(30 (10 + t)) + E^(20 (20 + t))
- 2 E^(10 (30 + t)) - 2 E^(100 + 30 t)) c[t]^2 x[t] == 330183/6250 + (53 x''[t])/25 *)
ic14x = Collect[ic14 /. sy, x, Simplify]
(* -(((1 + E^200) c x^2)/(10 E^200)) == 1/20 *)


Use the procedure described here to solve for x and c.

sn = NDSolveValue[{ode2x /.c -> c[t], ic14x /. c -> c[L], ic21, ic22, c'[t] == 0},
{x[t], c}, t, Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {x == -1/2, c == I/8}}] // Flatten;
sn//Last
(* 0.000769293 + 0.12592 I *)


which is c.

Plot[Evaluate@ReIm@First@sn, {t, 0, L},
ImageSize -> Large, AxesLabel -> {t, x}, LabelStyle -> {15, Bold, Black}]
Plot[Evaluate@ReIm@Last[y /. sy /. c -> Last[sn]], {t, 0, L},
ImageSize -> Large, AxesLabel -> {t, y}, LabelStyle -> {15, Bold, Black}]  Because two of the three conditions, {ic14x, ic21, ic22} are evaluated at x = L and only one at x = 0, integrating from L to 0 is both faster and more robust that integrating from 0 to L. This permits the rapid computation of four distinct eigenfunction-like families of solutions. To begin, solve ic14x for x[L] in terms of c[L].

sxL = (Solve[ic14x, x[L]] // Flatten) /. Rule -> Equal
(* {x == -((I E^100)/(Sqrt[2 + 2 E^200] Sqrt[c])),
x == (I E^100)/(Sqrt[2 + 2 E^200] Sqrt[c])} *)


Now, find and plot the first twelve solutions for sxL//First and Im[c] > 0.

Sort[Union[Flatten@Table[Quiet@Check[
NDSolveValue[{ode2x, sxL // First, ic21, ic22, c'[t] == 0}, {c[L]}, t,
Method -> {"Shooting", "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {c[L] == n I}}] // Flatten,
Nothing], {n, .105, .135, .0005}],
SameTest -> (Abs[#1 - #2] < 10^-4 &)], Im[#1] < Im[#2] &];
GraphicsGrid[Partition[(s =
NDSolveValue[{ode2x, sxL // First, ic21, ic22, c'[t] == 0}, {x[t], c[L]}, t,
Method -> {"Shooting", "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {c[L] == #}}];
Plot[Evaluate@ReIm@First@s, {t, 0, L}, PlotLabel -> Last@s]) & /@ %,
UpTo], ImageSize -> Large] The identical code with sxL // First replaced by sxL // Last yields The remaining two sets are obtained as above but with the "StartingInitialConditions" c[L] == n I replaced by c[L] == -n I.  I do not believe that there are other solution families.

To illustrate the accuracy of the computations above, first apply sy, to {ode1, ic11, ic12, ic13}

Simplify[{ode1, ic11, ic12, ic13} /. sy]
(* {True, True, True, True} *)


A typical evaluation of x, obtained using the procedure introduce in the preceding addendum,

sn = NDSolveValue[{ode2x, sxL // First, ic21, ic22, c'[t] == 0}, {x, c[L]}, t,
Method -> {"Shooting", "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"},
"StartingInitialConditions" -> {c[L] == 126 10^-3 I}}, InterpolationOrder -> All]
// Flatten]


which yields, when plotted, the first plot in the answer. Now, compute the corresponding numerical error,

(Subtract @@ ode2x) /. x -> First[sn] /. c[t] -> Last[sn];
Plot[Norm@%%, {t, 0, L}, PlotRange -> {0, 10^-5}, ImageSize -> Large,
AxesLabel -> {t, err}, LabelStyle -> {15, Bold, Black}] The largest error, 0.00156 at t = L, is a bit large but nonetheless satisfactory. (Reducing the maximum error by using a larger WorkingPrecision makes no visible difference to a plot of the solution.) Application of sn to the remaining boundary conditions shows that they too are well satisfied.

(Subtract @@ First[sxL]) /. x -> First[sn] /. c -> Last[sn]
(* -6.66134*10^-16 - 2.22045*10^-16 I *)
(Subtract @@ ic21) /. x -> First[sn] /. c -> Last[sn]
(* 3.49637*10^-7 + 1.73382*10^-8 I *)
(Subtract @@ ic22) /. x -> First[sn] /. c -> Last[sn]
(* 0. + 0. I *)

• Thanks ! As mentioned, the initial value of x needs to equal one. – ABCDEMMM Aug 13 at 3:02
• @ABCDEMMM You cannot force x == 1` unless you use it as a boundary condition, which you have not done. – bbgodfrey Aug 13 at 3:23
• to check the results, firstly, we can set ic14: x[L]*x[L]*(y'[L] - 0.01*y'''[L]) == 0.; And the correct results must be that: x = 1 and y = 0 for t =[0 10], however, the proposed method does not provid this result. – ABCDEMMM Aug 13 at 23:30
• we just want to check your method in MMA, then we use ic14: x[L]*x[L]*(y'[L] - 0.01*y'''[L]) == 0. if we use this new bc, then y must be zero, Next step, if y equals 0, x will never be changed in ode2, which means the resulting x must equal the initial value x_0. – ABCDEMMM Aug 14 at 0:48
• As I mentioned, for the test job (a simple test procedure), we use the new ic4: x[L]*x[L]*(y'[L] - 10^-2*y'''[L]) == 0, not old one 1/20. The verify process is that use this new bc, then y must be zero, Next step, if y equals 0, x will never be changed in ode2, which means the resulting x must equal the initial value x_0 – ABCDEMMM Aug 14 at 1:31