Solve Falkner Skan Numerically?

Question

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!

Answer 1

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);

end

And an implementation of the shooting method is:

clear

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);
    plot(x_out,y_out)
    pause

    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;
end

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.