[Date Index] [Thread Index] [Author Index]

Re: How to solve the Blasius eq. in Mathematica?

• To: mathgroup at smc.vnet.net
• Subject: [mg14122] Re: How to solve the Blasius eq. in Mathematica?
• From: lawry at maths.ox.ac.uk (James Lawry)
• Date: Fri, 25 Sep 1998 03:15:46 -0400
• Organization: Oxford Centre for Industrial & Applied Mathematics
• References: <6ud20f$16g@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Jason Krywicki <krywicki_NOSPAM_ at bucknell.edu> wrote:
>I was wondering what the best way would be to go about solving 
>
>   2f''' + ff'' == 0  f[0]==0, f'[0]==0,f'[inf]==1
>
>in Mathematica would be? I guess NDSolve, but is there anything special
>i have to do?

Yes. First, because the equation is nonlinear and the boundary
conditions are not all imposed at one point, the built-in NDSolve
cannot do the whole problem for you and you will need to use something
like a shooting method using NDSolve in combination with FindRoot:
effectively you guess a value of f''[0], solve the differential
equation with the three conditions at 0, and then use FindRoot to vary
the value of f''[0] until it matches the boundary condition at the
other end.

Second, the boundary condition imposed at infinity causes difficulty to
Mathematica. Two ways around this are: either transform the independent
variable (e.g. with t = 1 - Exp[-x]) so that x=infinity corresponds to
a finite value in the new variable; or else (the quick and dirty way)
decide that some finite value is "far enough away" to be considered
infinity. In the case of your equation, x=10 is amply far enough away
for several digits precision.

A simple implementation of a shooting routine is available from my
webpage at http://users.ox.ac.uk/~chri0220/software/DShoot.m

Using it, the command to solve your equation is
  soln = DShoot[{2 f'''[x] + f[x] f''[x] == 0, f[0] == 0, f'[0] == 0}, 
                f'[10] == 1, f, {x, 0, 10}, {f''[0], 0.3, 0.5}]

(where I have used 10 in place of infinity and given the range (0.3,
0.5) as two initial values for FindRoot to use to find the correct
value of f''[0]. The solution returned is 

  {{f -> InterpolatingFunction[{{0., 10.}}, <>]}, f''[0] == 0.332058}

so the root for f''[0] satisfying the condition f'[10] == 1 has been
found, and plotting f'[x] with

  Plot[Evaluate[(f /. soln[[1]])'[x]], {x, 0, 10}, PlotRange  -> {0, 1}]

shows the familar viscous boundary layer shape. (Always worth checking,
as using shooting and a finite boundary condition could produce
spurious results.)

James Lawry.