# Solving a boundary value problem

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


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.

• WOW!. I tried a very large number of initial guesses but no luck. Thanks for such a nice answer. – zhk Jan 13 '17 at 15:13
• Thanks, @MMM. I also tried many initial guesses without success until I used the alternative boundary conditions. – bbgodfrey Jan 13 '17 at 15:23
• Thank you so much for helping me to solve this kind of question! :) I am very new to this, hence i can only manage to solve simple question, not this kind of question. However, may I know how do you plot G[n] and H[n] then? We cant simply change f[n] = First@ans to g[n] = First@ans; right? @bbgodfrey – dayana nisa Feb 3 '17 at 1:07