2
$\begingroup$

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!

  • 0
    Woops! beta is a real number if it helps!2012-06-15
  • 2
    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 (http://en.wikipedia.org/wiki/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.2012-06-15
  • 0
    @in_wolfram_we_trust, sounds like an answer, after a bit more elaboration.2012-06-15
  • 0
    I'd elaborate, but I'm lazy. It'll do @Joe some good to work through the details himself.2012-06-15
  • 0
    @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!2012-06-15
  • 0
    @Joe, here's the solution I found. http://imgur.com/tbdac2012-06-15
  • 0
    @in_wolfram we trust Thank you so much for your thorough reply! I really appreciate it! I did something similar but your code is much more straightforward and structured!2012-07-10

1 Answers 1

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.

  • 0
    Thank you so much for the reply in_wolfram_we_trust!2012-07-10