Solving stiff equations in Mathematica

Question

I have problem to solve stiff equations. Any idea on how to solve this? I have tried "StiffSwitching" but it didnt work. Im solving this using Mathematica 10.

Here is my code. Im sorry if I wrote the code badly but I have fllowed the instructions on how to copy and paste code from Mathematica and this is what I end up with.

With[{n = 0.6}, ODE3[L_] := {Derivative[1][H][\[Eta]] == -2*F[\[Eta]] - ((1 - n)/(1 + n))*\[Eta]*Derivative[1][F][\[Eta]], 
        F[\[Eta]]^2 - (G[\[Eta]] + 1)^2 + (H[\[Eta]] + ((1 - n)/(1 + n))*\[Eta]*F[\[Eta]])*Derivative[1][F][\[Eta]] == 
         (Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 1)/2)*Derivative[2][F][\[Eta]] + 
          Derivative[1][F][\[Eta]]*((n - 1)*(Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 3)/2)*(Derivative[1][F][\[Eta]]*Derivative[2][F][\[Eta]] + 
             Derivative[1][G][\[Eta]]*Derivative[2][G][\[Eta]])), 2*F[\[Eta]]*(G[\[Eta]] + 1) + (H[\[Eta]] + ((1 - n)/(1 + n))*\[Eta]*F[\[Eta]])*Derivative[1][G][\[Eta]] == 
         (Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 1)/2)*Derivative[2][G][\[Eta]] + 
          Derivative[1][G][\[Eta]]*((n - 1)*(Derivative[1][F][\[Eta]]^2 + Derivative[1][G][\[Eta]]^2)^((n - 3)/2)*(Derivative[1][F][\[Eta]]*Derivative[2][F][\[Eta]] + 
             Derivative[1][G][\[Eta]]*Derivative[2][G][\[Eta]])), F[0] == 0, G[0] == 0, H[0] == 0, F[L] == 0, G[L] == 0}]


sol4[L_] := Flatten[NDSolve[ODE3[L], {F, G, H}, \[Eta], 
    Method -> {"Shooting", "StartingInitialConditions" -> {G[0] == 0, Derivative[1][G][0] == -0.65, F[0] == 0, Derivative[1][F][0] == 0.55, 
        H[0] == 0}}]]


Plot[Evaluate[{F[\[Eta]], G[\[Eta]], H[\[Eta]]} /. sol4[14]], {\[Eta],
   0, 14}, PlotStyle -> Thick, Frame -> True, 
 FrameLabel -> {"\[Eta]", "F, G, H"}, 
 Epilog -> {Dashed, Thickness[0.005], 
   Line[{{0, 0.8844}, {8, 0.8844}}]}]

I got this type of error:

At \[Eta] == 11.976704516552642`, step size is effectively zero; singularity or stiff system suspected. 

This only works if I let n=1. But I have to plot for different values of n ranging from 0.6 to 1.4. Below is what plots I got from Mathematica.

It would be great if someone could help me and explain to me. Thank you in advance!

Answer 1

When solving a boundary value problem, the shooting method uses nonlinear rootfinding to find the initial conditions that would satisfy the boundary conditions. To see what's going on here, replace your BVP with the IVP corresponding to the first guess that you told the shooting method to start with. That is, replace

F[L] == 0, G[L] == 0

with

F'[0] == 0.55, G'[0] == -0.65

and plot the resulting solutions of the IVP. Unlike when solving the BVP, Mathematica will actually give you the solution up to time 11.977 when the singularity occurs. You can clearly see in the plot how the functions blow up — so it is definitely not stiffness causing this:

This makes it clear that your problem is that certain (many?) initial values result in singularities. So the shooting method is not suitable, unless by some coincidence the initial values that it tries just happen to not lead to singularities.

When I tried setting the starting initial conditions to F'[0] == 0.5, G'[0] == -0.7 (a bit of guesswork and trial and error), it seemed to work correctly, giving a solution that satisfies the BCs: