This example shows how to use continuation to gradually extend a BVP solution to larger intervals.

Falkner-Skan boundary value problems [1] arise from similarity solutions of a viscous, incompressible, laminar flow over a flat plate. An example equation is

${{{\mathit{f}}^{\prime}}^{\prime}}^{\prime}+\mathit{f}\text{\hspace{0.17em}}{{\mathit{f}}^{\prime}}^{\prime}+\beta \left(1-{{\mathit{f}}^{\prime}}^{\text{\hspace{0.17em}}2}\right)=0$.

The problem is posed on the infinite interval $\left[0,\infty \right]$ with $\beta =0.5$, subject to the boundary conditions

$\mathit{f}\left(0\right)=0$,

${\mathit{f}}^{\prime}\left(0\right)=0$,

${\mathit{f}}^{\prime}\left(\infty \right)=1$.

The BVP cannot be solved on the infinite interval, and it is impractical to solve the BVP in a very large finite interval. Instead, this example solves a sequence of problems posed on the smaller interval $\left[0,\mathit{a}\right]$ to verify that the solution has consistent behavior as $\mathit{a}\to \infty $. This practice of breaking the problem up into simpler problems, with the solution of each problem feeding back in as the initial guess for the next problem, is called *continuation*.

To solve this system of equations in MATLAB, you need to code the equations, boundary conditions, and options before calling the boundary value problem solver `bvp4c`

. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Create a function to code the equations. This function should have the signature `dfdx = fsode(x,f)`

, where:

x is the independent variable.

f is the dependent variable.

You can rewrite the third-order equation as a system of first-order equations using the substitutions ${\mathit{f}}_{1}=\mathit{f}$, ${\mathit{f}}_{2}={\mathit{f}}^{\prime}$, and ${\mathit{f}}_{3}={{\mathit{f}}^{\prime}}^{\prime}$. The equations become

${{\mathit{f}}_{1}}^{\prime}={\mathit{f}}_{2}$,

${{\mathit{f}}_{2}}^{\prime}={\mathit{f}}_{3}$,

${{\mathit{f}}_{3}}^{\prime}=-{\mathit{f}}_{1}{\mathit{f}}_{3}-\beta \left(1-{\mathit{f}}_{2}^{2}\right)$.

The corresponding function is

function dfdeta = fsode(x,f) b = 0.5; dfdeta = [ f(2) f(3) -f(1)*f(3) - b*(1 - f(2)^2) ]; end

*Note: All functions are included as local functions at the end of the example.*

Now, write a function that returns the residual value of the boundary conditions at the boundary points. This function should have the signature `res = fsbc(f0,finf)`

, where:

`f0`

is the value of the boundary condition at the beginning of the interval.`finf`

is the value of the boundary condition at the end of the interval.

To calculate the residual values, you need to put the boundary conditions in the form $\mathit{g}\left(\mathit{x},\mathit{y}\right)=0$. In this form, the boundary conditions are

$\mathit{f}\left(0\right)=0$,

${\mathit{f}}^{\prime}\left(0\right)=0$,

${\mathit{f}}^{\prime}\left(\infty \right)-1=0$.

The corresponding function is

function res = fsbc(f0,finf) res = [f0(1) f0(2) finf(2) - 1]; end

Lastly, you must provide an initial guess for the solution. A crude mesh of five points and a constant guess that satisfies the boundary conditions are good enough to get convergence on the interval $\left[0,3\right]$. The variable `infinity`

denotes the right-hand limit of the interval of integration. As the value of `infinity`

increases on subsequent iterations from 3 to its maximum value of 6, the solution from each previous iteration acts as the initial guess for the next iteration.

infinity = 3; maxinfinity = 6; solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);

Solve the problem in the initial interval $\left[0,\text{\hspace{0.17em}}3\right]$. Plot the solution, and compare the value of ${{\mathit{f}}^{\prime}}^{\prime}\left(0\right)$ to the analytic value [1].

sol = bvp4c(@fsode,@fsbc,solinit); x = sol.x; f = sol.y; plot(x,f(2,:),x(end),f(2,end),'o'); axis([0 maxinfinity 0 1.4]); title('Falkner-Skan Equation, Positive Wall Shear, \beta = 0.5.') xlabel('x') ylabel('df/dx') hold on

`fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n')`

Cebeci & Keller report that f''(0) = 0.92768.

fprintf('Value computed using infinity = %g is %7.5f.\n', ... infinity,f(3,1))

Value computed using infinity = 3 is 0.92915.

Now, solve the problem on progressively larger intervals by increasing the value of `infinity`

for each iteration. The `bvpinit`

function extrapolates each solution to the new interval to act as the initial guess for the next value of `infinity`

. Each iteration prints the calculated value of ${{\mathit{f}}^{\prime}}^{\prime}\left(0\right)$ and superimposes the plot of the solution over the previous solutions. When `infinity = 6`

, the consistent behavior of the solution becomes evident and the value of ${{\mathit{f}}^{\prime}}^{\prime}\left(0\right)$ is very close to the predicted value.

for Bnew = infinity+1:maxinfinity solinit = bvpinit(sol,[0 Bnew]); % Extend solution to new interval sol = bvp4c(@fsode,@fsbc,solinit); x = sol.x; f = sol.y; fprintf('Value computed using infinity = %g is %7.5f.\n', ... Bnew,f(3,1)) plot(x,f(2,:),x(end),f(2,end),'o'); drawnow end

Value computed using infinity = 4 is 0.92774.

Value computed using infinity = 5 is 0.92770.

Value computed using infinity = 6 is 0.92770.

`hold off`

Listed here are the local helper functions that the BVP solver `bvp4c`

calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function dfdeta = fsode(x,f) % equation being solved dfdeta = [ f(2) f(3) -f(1)*f(3) - 0.5*(1 - f(2)^2) ]; end %------------------------------------------- function res = fsbc(f0,finf) % boundary conditions res = [f0(1) f0(2) finf(2) - 1]; end %-------------------------------------------

[1] Cebeci, T. and H. B. Keller. "Shooting and Parallel Shooting Methods for Solving the Falkner-Skan Boundary-layer Equation*.*" *J. Comp. Phys.*, Vol. 7, 1971, pp. 289-300.