Solving a boundary value problem

Question

I have a problem solving a boundary layer problem with an infinity boundary conditions. I am almost there (I think). The problem is that, mathematica gives me this error:

NDSolve: At η == 6.53483 , step size is effectively zero; singularity or stiff system suspected

The system is exactly what I need to solve, but maybe my shooting vaue is not correct or something else is wrong.

This is the graph that I should get. The graph is F[η] over η for different n:

I only try n = 1.0 first, but still the shape of the graph is not even same. Can anyone help me?

This is my code:

inf = 20; 
sol = 
  NDSolve[
    {H'[η] == -2*F[η] - ((1 - n)/(1 + n))*η*F'[η], 
     F[η]^2 - (G[η] + 1)^2 + (H[η] + ((1 - n)/(1 + n))*η*F[η])*F'[η] ==  
       (F'[η]^2 + G'[η]^2)^((n - 1)/2)*F''[η] + 
       F'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)*(F'[η]*F''[η] + G'[η]*G''[η])), 
     2*F[η]*(G[η] + 1) + (H[η] + ((1 - n)/(1 + n))*η*F[η])*G'[η] == 
       (F'[η]^2 + G'[η]^2)^((n - 1)/2)*G''[η] + 
       G'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)*(F'[η]*F''[η] + G'[η]*G''[η])),  
     F[0] == 0, G[0] == 0, H[0] == 0, F[inf] == 0, G[inf] ==- 1} /. n -> 1.0,
    {F,G,H}, {η, 0, inf}, 
    Method -> 
     {"Shooting", "StartingInitialConditions" -> { F'[0]== 0.5,G'[0] ==-0.7}}];

Plot[Evaluate[{F[η]} /. sol], {η,0.00001,inf}, 
  PlotRange -> All, AxesLabel -> {x, y}]

Answer 1

Solution for n == 1

As noted by xzczd, this problem is similar to that in 104170, although more complicated. The answer is easy to obtain for n == 1, if we happen to know good initial guesses for F'[0] and G'[0]. For instance, using {F'[0] == .510, G'[0] == -0.616} immediately gives

The complete solution is

Plot[Evaluate[{F[η],G[η],H[η]} /. sol], {η, 0, inf}, PlotRange -> All, AxesLabel -> {x, y}]

The question now becomes, how to obtain these good initial guesses. Proceed as follows. First, explicitly set n == 1 in the second and third ODEs, solve for the second derivatives, and Simplify.

(Flatten@Solve[{F[η]^2 - (G[η] + 1)^2 + (H[η] + ((1 - n)/(1 + n))*η*F[η])*F'[η] ==  
    (F'[η]^2 + G'[η]^2)^((n - 1)/2)*F''[η] + 
    F'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)*(F'[η]*F''[η] + G'[η]*G''[η])), 
    2*F[η]*(G[η] + 1) + (H[η] + ((1 - n)/(1 + n))*η*F[η])*G'[η] == 
    (F'[η]^2 + G'[η]^2)^((n - 1)/2)*G''[η] + 
    G'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)*(F'[η]*F''[η] + G'[η]*G''[η]))}] 
    // Simplify) /. Rule -> Equal
(* {F''[η] == -1 + F[η]^2 - 2 G[η] - G[η]^2 + H[η] F'[η],
    G''[η] == 2 F[η] (1 + G[η]) + H[η] G'[η]} *)

Also, the first ODE reduces to

H'[η] == -2*F[η]

Next, derive alternative boundary conditions at large η. This can be done by linearizing the equations about F[η] == 0 and 1 + G[η] == 0,

{F''[η] == H[η] F'[η], G''[η] == H[η] G'[η]} 

Now, it is evident from the first ODE that H is approximately constant for large η, allowing the linearized equations to be integrated once (with constants of integration chosen so that F'[η] and G'[η] approach zero at large η.

{F'[η] == H[η] F[η], G'[η] == H[η] (G[η] + 1)} 

The critical advantage of these new boundary conditions at inf is that they are valid for fairly small values of inf. Combine these n == 1 ODEs with the new boundary conditions and choose inf == 5 to start.

inf = 5; 
sol = NDSolve[{H'[η] == -2*F[η] , 
    F''[η] == -1 + F[η]^2 - 2 G[η] - G[η]^2 + H[η] F'[η],
    G''[η] == 2 F[η] (1 + G[η]) + H[η] G'[η], 
    F[0] == 0, G[0] == 0, H[0] == 0, 
    F'[inf] == H[inf] F[inf], G'[inf] == H[inf] (G[inf] + 1)}, 
    {F, G, H}, {η, 0, inf}, Method -> {"Shooting", "StartingInitialConditions" 
    -> {F[0] == 0, G[0] == 0, H[0] == 0, F'[0] == .5, G'[0] == -0.7}}];
Flatten@{F'[0], G'[0]} /. sol
(* {0.51028, -0.615914} *)

where the initial guess for {F'[0], G'[0]} are those in the question. The result, immediately above, is a better set of initial guesses. So, use them with inf == 20, yielding still better initial guesses (although the ones just above are good enough).

(* {0.510233, -0.615922} *)

which have hardly changed. These are the values on which the initial guesses for the computation at the beginning of the answer are based.

Solutions for various n (corrected)

Solutions for other values of n can be obtained iteratively by using the initial guesses for n == 1 for a slightly smaller value of n, and so forth. However, inf must be somewhat larger to assure that H[η] is indeed constant near η == inf Start by defining

inf = 50;
s = ParametricNDSolve[{H'[η] == -2*F[η] - ((1 - n)/(1 + n))*η*F'[η], 
    F[η]^2 - (G[η] + 1)^2 + (H[η] + ((1 - n)/(1 + n))*η*F[η])*F'[η] == 
    (F'[η]^2 + G'[η]^2)^((n - 1)/2)*F''[η] + F'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)
    *(F'[η]*F''[η] + G'[η]*G''[η])), 
    2*F[η]*(G[η] + 1) + (H[η] + ((1 - n)/(1 + n))*η*F[η])*G'[η] == 
    (F'[η]^2 + G'[η]^2)^((n - 1)/2)*G''[η] + G'[η]*((n - 1)*(F'[η]^2 + G'[η]^2)^((n - 3)/2)
    *(F'[η]*F''[η] + G'[η]*G''[η])), 
    F[0] == 0, G[0] == 0, H[0] == 0, 
    F'[inf] == H[inf] F[inf], G'[inf] == H[inf] (G[inf] + 1)}, 
    {F, G, H, F', G'}, {η, 0, inf}, {n, Fp0, Gp0}, 
    Method -> {"Shooting", "StartingInitialConditions" -> 
    {F[0] == 0, G[0] == 0, H[0] == 0, F'[0] == Fp0, G'[0] == Gp0}}];

Then, slowly decrease n, using the computed values of {F'[0], G'[0]} from a value of n as initial guesses for the next value of n. It turns out that step in n must be quite small as n approaches 0.6.

guess = {0.510233, -0.615922};
Do[ans = Through[{F, G, H}[n, First@guess, Last@guess]] /. s;
    f[n] = First@ans;
    guess = D[Through[ans[η]], η][[1 ;; 2]] /. η -> 0;,
    {n, 1., .6, -.002}]
Plot[Evaluate@Through[Thread[f[n] /. n -> Range[1., .6, -.1]][η]], {η, 0, 20}, 
    PlotRange -> All, AxesLabel -> {x, y}, AspectRatio -> 5/7]

The total computation requires about 85 sec on my PC.