5
$\begingroup$

Here are the code that I had tried but there is an error

NDSolve[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0,    
    f'[0] == 0, f'[y -> ∞] -> 2}, f, {y, 0, 10}]

Here are the question :

enter image description here

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Commented Nov 28, 2015 at 2:51
  • 1
    $\begingroup$ What is the symbol Pr? Note also that NDSolve cannot operate over an infinite domain. So, the upper limit on y must be some finite number. $\endgroup$
    – bbgodfrey
    Commented Nov 28, 2015 at 5:00
  • 1
    $\begingroup$ I have a feeling this kind of equation has been asked several times in this site. Can you add some background information of the equation so we might be able to solve the problem better? $\endgroup$
    – xzczd ♦
    Commented Nov 28, 2015 at 6:39
  • 1
    $\begingroup$ @xzczd Certainly. I shall do so tomorrow morning. Just now, I would observe that such problems typically are very difficult to solve numerically, because the solution tracks an unstable attractor of the nonlinear equation. This one is relatively easy, however, because it enters the asymptotic regime so quickly. $\endgroup$
    – bbgodfrey
    Commented Nov 28, 2015 at 7:06

2 Answers 2

Reset to default
3
$\begingroup$

f can be computed as follows.

s = NDSolveValue[{f'''[y] + f[y] f''[y] + 4 - (f'[y])^2 == 0, f[0] == 0, f'[0] == 0, 
    f'[5] == 2}, f, {y, 0, 5}, Method -> "Shooting", 
    "StartingInitialConditions" -> {f[0] == 0, f'[0] == 0, f''[0] == 3.48}}];
Plot[s[y], {y, 0, 5}, AxesLabel -> {y, f}]
Plot[s'[y], {y, 0, 5}, PlotRange -> {-1, 3}, AxesLabel -> {y, f'}]

enter image description here enter image description here

Clearly, the solution already is in the asymptotic regime by y = 3.

Computation of θ

With ψ defined as θ', the second ODE can be solved formally as 

DSolve[{ψ'[y]/ψ[y] + Pr f[y] == 0, ψ[0] == -1}, ψ[y], y]
{{ψ[y] -> -E^(-Integrate[-(Pr*f[K[1]]), {K[1], 1, 0}] + 
           Integrate[-(Pr*f[K[1]]), {K[1], 1, y}])}}

To proceed, Pr needs to be specified. For now, set it to unity. Then, the answer just obtained can be written as the following function.

t[y_?NumericQ] := -Exp[- NIntegrate[s[yp], {yp, 0, y}, 
      Method -> {Automatic, "SymbolicProcessing" -> False}]];
Plot[t[y], {y, 0, 5}, AxesLabel -> {y, θ'}]

enter image description here

Finally, θ can be obtained by a second integration, with c5 the constant of integration chosen so that θ[5] == 0, as required by the question.

c5 = NIntegrate[t[yp], {yp, 0, 5}]
(* -1.23953 *)
tt[y_] := NIntegrate[t[yp], {yp, 0, y}] - c5;
Plot[tt[y], {y, 0, 5}, AxesLabel -> {y, θ}, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ How you chose, 'f''[0] == 3.48'? $\endgroup$
    – zhk
    Commented Nov 28, 2015 at 6:38
  • 1
    $\begingroup$ @MMM 3.48 is only a guess, so I does not need to be precise. I tried a few guesses until one worked. $\endgroup$
    – bbgodfrey
    Commented Nov 28, 2015 at 6:52
1
$\begingroup$

@bbgodfrey has explained very well his approach in the answer to this question. But, In my answer, I am trying to solve both the equations simultaneously. 

Eqn1 = f'''[x] + f[x] f''[x] + 4 - (f'[x])^2 == 0

Eqn2 = T''[x] + Pr f[x] T'[x] == 0

BC1 = f[0] == 0;

BC2 =  f'[0] == 0;

BC3 = f'[inf1] == 2;

BC4 = T'[0] == -1;

BC5 = T[inf1] == 0;

param1 = {Pr -> 3.97};

inf1 = 5;

Sol1 = NDSolve[{Eqn1, Eqn2, BC1, BC2, BC3, BC4, BC5} /. param1, {f, 
   T}, {x, 0, inf1}, Method -> {"Shooting", 
"StartingInitialConditions" -> {f[0] == 0, f'[0] == 0, f''[0] == 3.48,
 T[0] == 0, T'[0] == -10}}]

Plot[{f'[x] /. Sol1, T[x] /. Sol1}, {x, 0, inf1}, PlotRange -> All, 
 PlotStyle -> {Black, Red}, Frame -> True, 
 FrameStyle -> Directive[Black, Bold, 12], PlotRange -> All, 
 Axes -> False]

enter image description here

I hope this helps.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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