3
$\begingroup$

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] == 0;
ic12 = y''[0] == 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] == 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[0] == 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] == 0;
ic12 = y''[0] == 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] == 0;
ic22 = x'[L] == 0;
sy = (DSolve[{ode1, ic11, ic12, ic13}, y, t] // Flatten) /. 
  C[1] -> c Exp[-100]
ode2x = Simplify[ode2 /. sy];

ic14x = Collect[ic14 /. sy, x[10], Simplify];
sn = NDSolveValue[{ode2x /. c -> c[t], ic14x /. c -> c[L], ic21, ic22,
      c'[t] == 0}, {x[t], c[10]}, t, 
    Method -> {"Shooting", 
      "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}, 
      "StartingInitialConditions" -> {x[0] == -1/2, c[0] == 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:

enter image description here

Obviously, x is not const, oscillation can be observed.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Aug 12, 2020 at 22:21
  • $\begingroup$ 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.] $\endgroup$
    – Michael E2
    Commented Aug 12, 2020 at 22:23
  • $\begingroup$ @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. $\endgroup$
    – ABCDEMMM
    Commented Aug 12, 2020 at 23:04
  • $\begingroup$ @ABCDEMMM Requiring that x[0] == 1 as a boundary condition would overspecify this ode system. Note, however, that "StartingInitialConditions" -> {x[0] == 1} does not guarantee that the solution for x is equal to 1 at t = 0. $\endgroup$
    – bbgodfrey
    Commented Aug 12, 2020 at 23:59

1 Answer 1

Reset to default
5
$\begingroup$

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[1] -> 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[10], Simplify]
(* -(((1 + E^200) c x[10]^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[10]}, t, Method -> {"Shooting", 
    "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}, 
    "StartingInitialConditions" -> {x[0] == -1/2, c[0] == 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}]

enter image description here

enter image description here

Addendum: Four families of solutions

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[10] == -((I E^100)/(Sqrt[2 + 2 E^200] Sqrt[c[10]])), 
    x[10] == (I E^100)/(Sqrt[2 + 2 E^200] Sqrt[c[10]])} *)

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[3]], ImageSize -> Large]

enter image description here

The identical code with sxL // First replaced by sxL // Last yields

enter image description here

The remaining two sets are obtained as above but with the "StartingInitialConditions" c[L] == n I replaced by c[L] == -n I.

enter image description here

enter image description here

I do not believe that there are other solution families.

Second Addendum: Accuracy Demonstration

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

enter image description here

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[10] -> Last[sn]
(* -6.66134*10^-16 - 2.22045*10^-16 I *)
(Subtract @@ ic21) /. x -> First[sn] /. c[0] -> Last[sn]
(* 3.49637*10^-7 + 1.73382*10^-8 I *)
(Subtract @@ ic22) /. x -> First[sn] /. c[10] -> Last[sn]
(* 0. + 0. I *)
Improve this answer
$\endgroup$
13
  • $\begingroup$ Thanks ! As mentioned, the initial value of x needs to equal one. $\endgroup$
    – ABCDEMMM
    Commented Aug 13, 2020 at 3:02
  • $\begingroup$ @ABCDEMMM You cannot force x[0] == 1 unless you use it as a boundary condition, which you have not done. $\endgroup$
    – bbgodfrey
    Commented Aug 13, 2020 at 3:23
  • $\begingroup$ 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. $\endgroup$
    – ABCDEMMM
    Commented Aug 13, 2020 at 23:30
  • $\begingroup$ 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. $\endgroup$
    – ABCDEMMM
    Commented Aug 14, 2020 at 0:48
  • $\begingroup$ 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 $\endgroup$
    – ABCDEMMM
    Commented Aug 14, 2020 at 1:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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