Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

how would we go about solving the Falkner Skan equation numerically?

The equation is

$$f'''+ff''+\beta\left(1-f'^2\right)=0$$ $$f(0) = f'(0) = 0$$ $$f'(\infty) = 1$$

I tried using ODE45 to solve this but realised an $f''(0)$ condition is lacking and moreover, how is one supposed to deal with $f'(\infty) = 1$?

Thanks in advance for any incoming comments!

share|cite|improve this question
Woops! beta is a real number if it helps! – Joe Jun 15 '12 at 10:58
ODE45 implements a Runge-Kutta method which solves initial value problems. As you've realised this is not an initial value problem, it's a boundary value problem on the half-line. The simplest approach is probably the shooting method (, which reduces a BVP to an IVP. Obviously you can't integrate to infinity numerically so you'll have to pick a sufficiently big $L$ and treat that as your right hand endpoint. – in_wolframAlpha_we_trust Jun 15 '12 at 11:12
@in_wolfram_we_trust, sounds like an answer, after a bit more elaboration. – J. M. Jun 15 '12 at 11:21
I'd elaborate, but I'm lazy. It'll do @Joe some good to work through the details himself. – in_wolframAlpha_we_trust Jun 15 '12 at 11:25
@in_wolfram_we_trust Yes! I only realised it a while ago that the problem is one that is a BVP instead of an IVP! I did suspect the use of Shooting Method is needed and your reply pretty much confirms it though I must say, I am not very familiar with the method nor am I familiar with coding it up in Matlab. Something for me to learn I guess! – Joe Jun 15 '12 at 11:30
up vote 1 down vote accepted

This is a boundary value problem (BVP), not an initial value problem (IVP). We can solve it using MATLAB's ODE45 solver by converting it from a BVP to an IVP through the shooting method. The general idea is the following:

  1. Disregard the boundary value(s) on the right. For this problem it means we ignore the condition that $f'(\infty) = 1.$
  2. Introduce new condition(s) on the left. In our case we've thrown out one condition on the right, so we have to introduce one condition on the left, say $f''(0) = 1$.
  3. Solve the initial value problem using the condition we introduced above. If we're very lucky the solution at the right hand endpoint agrees with the boundary conditions we discarded in step 1. More likely the solution will over- or undershoot the desired value.
  4. Make a sensible adjustment to the condition(s) we've introduced on the left and repeat from step 3. To decide what sort of adjustment is sensible you use your favourite root finding technique: bisection method, secant method, Newton's method etc.

Here are a few points (and some MATLAB code) for the implementation of the shooting method to this problem. Firstly, we need to rewrite our equation $$ f''' + ff'' + \beta(1-f'^2) =0$$ as a system of first order equations. This is not too difficult: $$ f' = g, $$ $$ g' = h, $$ $$ h' = -fh - \beta(1-g^2).$$ The corresponding function in MATLAB to be used with the ODE45 solver is something like:

function y_out = falkner_skan(x,y_in)

global beta

y_out(1,1) = y_in(2);
y_out(2,1) = y_in(3);
y_out(3,1) = -y_in(1)*y_in(3) - beta*(1-y_in(2)^2);


And an implementation of the shooting method is:


global beta
beta = 1;
right_hand_endpoint = 5;
stopping_distance = 0.01;

shoot_parameter(1) = 1;
initial_value = [0,0,shoot_parameter];
[x_out, y_out] = ode45(@falkner_skan,[0, right_hand_endpoint],initial_value);
f_dash_endpoint(1) = y_out(end,2);

shoot_parameter(2) = 1.1;

n = 2;
while abs(shoot_parameter(n) - shoot_parameter(n-1)) > stopping_distance 
    initial_value = [0,0,shoot_parameter(n)];
    [x_out, y_out] = ode45(@falkner_skan,[0, right_hand_endpoint],initial_value);
    f_dash_endpoint(n) = y_out(end,2);

    shoot_parameter(n+1) = shoot_parameter(n) - (f_dash_endpoint(n)-1)*(shoot_parameter(n)-shoot_parameter(n-1))/(f_dash_endpoint(n)-f_dash_endpoint(n-1));
    n = n+1;

In the code above I've used the secant method to make my sensible guesses and instead of solving from $0$ to infinity I've chosen a value of $5$ as my right hand endpoint. The code is uncommented, but it's only 25 lines long and there is nothing mysterious going on in it so you should be able to follow it fairly easily. Let me know if you have questions.

share|cite|improve this answer
Thank you so much for the reply in_wolfram_we_trust! – Joe Jul 10 '12 at 17:50

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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