4
$\begingroup$

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:

enter image description here

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}]
$\endgroup$
6
$\begingroup$

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

enter image description here

The complete solution is

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

enter image description here

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]

enter image description here

The total computation requires about 85 sec on my PC.

$\endgroup$
  • $\begingroup$ WOW!. I tried a very large number of initial guesses but no luck. Thanks for such a nice answer. $\endgroup$ – zhk Jan 13 '17 at 15:13
  • 1
    $\begingroup$ Thanks, @MMM. I also tried many initial guesses without success until I used the alternative boundary conditions. $\endgroup$ – bbgodfrey Jan 13 '17 at 15:23
  • $\begingroup$ 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 $\endgroup$ – dayana nisa Feb 3 '17 at 1:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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