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.