4 Runge-Kutta solution

In order to solve O.D.E’s such as the Blasius equation (5) we often need to resort to computer methods. Let us start by thinking about what an O.D.E actually represents. A first order O.D.E is a statement that the gradient of y, dy/dx, takes some value or function. We can write this as

dy-= f(x,y) dx

Given a particular starting point - some value of y corresponding to a given value of x - this equation specifies a particular curve in the x - y plane. A typical such solution is shown in figure 3.

Figure 3:

Typical 1st order O.D.E.

We could trace the curve defined by this O.D.E - i.e. solve the equation - by starting at this point and tracing forwards/backwards in steps of size h

e.g. dy-= y y(0) = 1 dx

In fact this is the basis for all numerical methods for solving O.D.E’s. Used like this, it is referred to as Euler’s method. Symbolically, we are writing the derivative as

dy-= yn+1---yn= f (x,y) dx h

which we can rewrite to find point y_n+1 from point y_n

yn+1 = yn + hf(xn,yn)

(6)

Obviously the smaller the step h being used, the more accurate the resulting solution will be. In fact, a more exact analysis shows that equation 6 is missing a term of the form h². We therefore say that the scheme is first order accurate. However using smaller steps h also increases the computational cost of the calculation, and it may be better to try to find more accurate schemes. It is possible to combine Euler steps of various sizes in a scheme in such a way as to produce a more accurate final scheme. The Runge-Kutta family of numerical schemes is constructed in this way. To solve the Blasius equation we will make use of the 4th order Runge-Kutta method, so called because it is 4th order accurate (the missing terms in the scheme are of the form h⁵). The 4th order Runge-Kutta method is available in MathCad in the form of a function call rkfixed. You will make use of this in Problem Sheet 3153-1 to solve the Blasius problem.

Of course, in this case we can solve the problem analytically. The solution is

x y = e

4.1 Higher order O.D.E’s

The Blasius equation is a 3rd order O.D.E, and the discussion above talked about 1st order O.D.E’s. It is possible to use Runge-Kutta for higher order O.D.E’s, but we must first manipulate the equation we want to solve. For instance, suppose we want to solve the equation

d2y (dy )2 --2 + --- = 0 dx dx

For this equation, define a new variable

dy- y1 = dx

Using this variable

d2y= dy1 dx2 dx

and we can write the 2nd order ODE as 2 1st order equations :

dy-= y1, dy1+ y21 = 0 dx dx

If we define a vector

( y) Y = y1

then we could write these as

dY ( y1) -dx = - y2 = f(Y,x) 1

which is in the same form as before - a 1st order O.D.E, but involving a vector of y values.

For the Blasius equation,

3 2 d-f + f-df-= 0 dz3 2 dz2

if we write

z = X, f = Y0, f '= Y1, f” = Y2

we solve the equations

dY0-= Y , dY1-= Y , dY2-= -1Y Y dX 1 dX 2 dX 2 0 2

Of course, we still need to find the boundary conditions in order to be able to solve these equations. We note that at the plate y = 0 the velocity is zero u_x = 0, whilst a long way from the plate the flow has to match the undisturbed flow u_x

U :

Boundary 1: At the plate y = 0, i.e. = 0, u_x = 0 i.e. f' = 0, f = 0.
Boundary 2: A long way from the plate, i.e. y , i.e. , u_x U, i.e. f' 1

This gives us a sufficient number of boundary conditions, but they are in two different locations. The rkfixed routine in MathCad requires you to specify values for Y ₀, Y ₁ and Y ₂ for a specific value of X, and will return the solution for other values of X starting at this point. The only way to find the solution in this case is to use a shooting method. We know that at X = 0, Y ₀ = Y ₁ = 0. Guess a value of Y ₂ at this point and then adjust it until the second boundary condition is satisfied.

[prev] [prev-tail] [front] [up]