# 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, G == 0, H == 0, F[inf] == 0, G[inf] ==- 1} /. n -> 1.0,
{F,G,H}, {η, 0, inf},
Method ->
{"Shooting", "StartingInitialConditions" -> { F'== 0.5,G' ==-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' and G'. For instance, using {F' == .510, G' == -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, G == 0, H == 0,
F'[inf] == H[inf] F[inf], G'[inf] == H[inf] (G[inf] + 1)},
{F, G, H}, {η, 0, inf}, Method -> {"Shooting", "StartingInitialConditions"
-> {F == 0, G == 0, H == 0, F' == .5, G' == -0.7}}];
Flatten@{F', G'} /. sol
(* {0.51028, -0.615914} *)


where the initial guess for {F', G'} 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, G == 0, H == 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, G == 0, H == 0, F' == Fp0, G' == Gp0}}];


Then, slowly decrease n, using the computed values of {F', G'} 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, 2017 at 15:13
• Thanks, @MMM. I also tried many initial guesses without success until I used the alternative boundary conditions. Jan 13, 2017 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 Feb 3, 2017 at 1:07