• Author / Uploaded
• malekan2005

Citation preview

5568ch20.qxd_jd

3/24/03

9:09 AM

Page 271

20 Solving Differential Equations

D

ifferential equations, those that define how the value of one variable changes with respect to another, are used to model a wide range of physical processes. You will use differential equations in chemistry, dynamics, fluid dynamics, thermodynamics, and almost every other scientific or engineering endeavor. A differential equation that has one independent variable is called an ordinary differential equation or ODE. Examples of ODEs include the equations to model the motion of a spring or the boundary layer equations from fluid dynamics. A partial differential equation (PDE) has more than one independent variable. The Navier-Stokes equations are an example of a set of coupled partial differential equations used in fluid dynamic analysis to represent the conservation of mass, momentum, and energy. This chapter will focus primarily on how to solve ordinary differential equations and will touch upon the more difficult to solve partial differential equations only briefly at the end of the chapter. We will discuss the difference between initial value and two-point boundary problems. We will write a class that represents a generic ODE and write two subclasses that represent the motion of a damped spring and a compressible boundary layer over a flat plate. We will develop a class named ODESolver that will define a number of methods used to solve ODEs and compare results generated by these methods with results from other sources. The specific topics covered in this chapter are—

271

5568ch20.qxd_jd

3/24/03

9:09 AM

272

Page 272

Chapter 20 Solving Differential Equations

• • • • • • • • • • • •

Ordinary differential equations The ODE class Initial value problems Runge-Kutta schemes Example problem: damped spring motion Embedded Runge-Kutta solvers Other ODE solution techniques Two-point boundary problems Shooting methods Example problem: compressible boundary layer Other two-point boundary solution techniques Partial differential equations

Ordinary Differential Equations An ODE is used to express the rate of change of one quantity with respect to another. You have probably been working with ODEs since you began your scientific or engineering course work. One defining characteristic of an ODE is that its derivatives are a function of one independent variable. A general form of a first-order ODE is shown in Eq. (20.1). dy + a1x2y + b1x2 + c = 0 dx

(20.1)

The order of a differential is defined as the order of the highest derivative appearing in the equation. Ordinary differential equations can be of any order. A general form of a second-order ODE is shown in Eq. (20.2). dy d2y + b1x2y + c1x2 + d = 0 2 + a1x2 dx dx

(20.2)

Any higher-order ODE can be expressed as a coupled set of first-order differential equations. For example, the second-order ODE shown in Eq. (20.2) can be reduced to a coupled set of two first-order differential equations. dy d dy a b = -a1x2 - b1x2y - c1x2 - d dx dx dx

(20.3)

dy d 1y2 = dx dx The second expression in Eq. (20.3) looks trivial in that the left-hand side is the same as the right-hand side, but the ODE solvers we will discuss later in

5568ch20.qxd_jd

3/24/03

9:09 AM

The ODE Class

Page 273

273

this chapter use the coupled first-order form of the ODE in their solution process. The ODE solvers would integrate the first-order equations shown in Eq. (20.3) to obtain values for the dependent variables y and dy/dx as a function of the independent variable x.

The ODE Class As you certainly know by this time, everything in Java is defined within a class. If we are working with ODEs we need to define a class that will encapsulate an ODE. We will write the ODE class to represent a generic ODE. It will be the superclass for specific ODE subclasses. The ODE class will declare fields and methods used by all ODE classes. Since an ODE is a mathematical entity, we will place the ODE class in the TechJava.MathLib package. When writing a class you must always consider the state and behavior of the item you are modeling. Let us first consider the fields that will define the state of an ODE. The ODE class will represent its associated ODE by one or more first-order differential equations. The ODE class will declare a field to store the number of first-order equations. Another field is needed to store the number of free variables in the ODE. Free variables are those that are not specified by boundary conditions at the beginning of the integration range. For initial value problems, the number of free variables will be zero. Two-point boundary problems will have one or more free variables. The coupled set of first-order ODEs is solved by integrating each of the ODEs step-wise over a certain range. The values of the independent and dependent variables will have to be stored at every step in the integration. To facilitate this, the ODE class declares two arrays— x[]which stores the values of the independent variable at each step of the integration domain and y[][] which stores the dependent variable or variables. The y[][] array is 2-D because an ODE might represent a system of first-order differential equations and therefore have more than one dependent variable. The ODE class constructor will take two input arguments that specify the number of first-order differential equations and number of free variables used by the ODE. Because the required number of steps along the integration path is not a fixed value, the x[] and y[][] arrays are allocated to a maximum number of steps. This approach may waste a little memory but is the simplest way to do things. Now let’s turn to the behavior of an ODE class. What does an ODE class have to do? It must declare a method to return the right-hand sides of the first-

5568ch20.qxd_jd

3/24/03

9:09 AM

274

Page 274

Chapter 20 Solving Differential Equations

order differential equations that describe the ODE. The ODE class will declare methods to return the number or first-order equations and free variables as well as methods to return the values of the x[] and y[][] arrays. There will be one method to return the entire array and another to return a single element of the array. The ODE class also declares methods to set the conditions at the start of the integration range and to compute the error at the end. These methods and the right-hand side method are ODE-specific. Since the ODE class represents a generic ODE, they are implemented as stubs. ODE subclasses will override these methods according to their needs. The ODE class code listing is shown next. package TechJava.MathLib; public class ODE { // This is used to allocate memory to the // x[] and y[][] arrays public static int MAX_STEPS = 999; // // // // //

numEqns = number of 1st order ODEs to be solved numFreeVariables = number of free variables at domain boundaries x[] = array of independent variables y[][] = array of dependent variables

private int numEqns, numFreeVariables; private double x[]; private double y[][]; public ODE(int numEqns, int numFreeVariables) { this.numEqns = numEqns; this.numFreeVariables = numFreeVariables; x = new double[MAX_STEPS]; y = new double[MAX_STEPS][numEqns]; } // //

These methods return the values of some of the fields.

public int getNumEqns() { return numEqns; } public int getNumFreeVariables() { return numFreeVariables; }

5568ch20.qxd_jd

3/24/03

9:09 AM

Page 275

Initial Value Problems

275

public double[] getX() { return x; } public double[][] getY() { return y; } public double getOneX(int step) { return x[step]; } public double getOneY(int step, int equation) { return y[step][equation]; } // //

This method lets you change one of the dependent or independent variables

public void setOneX(int step, double value) { x[step] = value; } public void setOneY(int step, int equation, double value) { y[step][equation] = value; } // //

These methods are implemented as stubs. Subclasses of ODE will override them.

public void getFunction(double x, double dy[], double ytmp[]) {} public void getError(double E[], double endY[]) {} public void setInitialConditions(double V[]) {} }

Initial Value Problems Before we discuss how to solve them, let’s explore a little bit about the nature of ODEs themselves. There are two basic types of boundary condition categories for ODEs—initial value problems and two-point boundary value problems. With an initial value problem, values for all of the dependent variables are specified at the beginning of the range of integration. The initial boundary serves as the “anchor” for the solution. The solution is marched outward from

5568ch20.qxd_jd

276

3/24/03

9:09 AM

Page 276

Chapter 20 Solving Differential Equations

the initial boundary by integrating the ODE at discrete steps of the independent variable. The dependent variables are computed at every step. Initial value problems are simpler to solve because you only have to integrate the ODE one time. The solution of a two-point boundary value problem usually involves iterating between the values at the beginning and end of the range of integration. The most commonly used techniques to solve initial value problem ODEs are called Runge-Kutta schemes and will be discussed in the next section.

Runge-Kutta Schemes One of the oldest and still most widely used groups of ODE integration algorithms is the Runge-Kutta family of methods. These are step-wise integration algorithms. Starting from an initial condition, the ODE is solved at discrete steps over the desired integration range. Runge-Kutta techniques are robust and will give good results as long as very high accuracy is not required. RungeKutta methods are not the fastest ODE solver techniques but their efficiency can be markedly improved if adaptive step sizing is used. Runge-Kutta methods are designed to solve first-order differential equations. They can be used on a single first-order ODE or on a coupled system of first-order ODEs. If a higher-order ODE can be expressed as a coupled system of first-order ODEs, Runge-Kutta methods can be used to solve it. To understand how Runge-Kutta methods work, consider a simple first-order differential equation. dy = f1x,y2 (20.4) dx To solve for the dependent variable y, Eq. (20.4) can be integrated in a step-wise manner. The derivative is replaced by its delta-form and the x term is moved to the right-hand side. ¢y = yn + 1 - yn = ¢xf1x,y2

(20.5)

Eq. (20.5) is the general form of the equation that is solved. Starting at an independent variable location xn where the value of yn is known, the value of the dependent variable at the next location, yn1, is equal to its value at the current location, yn, added to the independent variable step size, x, times the right-hand side function.

5568ch20.qxd_jd

3/24/03

9:09 AM

Page 277

277

Runge-Kutta Schemes

There is one question left to be resolved—where should we evaluate the right-hand side function? With the Euler method, as shown in Eq. (20.6), the function is evaluated at the current location, xn. yn + 1 = yn + ¢xf1xn,yn2

(20.6)

The value of y at the next step is computed using the slope of the f(x,y) function at the current step. If you perform a Taylor series expansion on Euler’s method you will find that it is first-order accurate in x. The Euler method is really only useful for linear or nearly linear functions. What happens, for instance, if the slope of the f(x,y) curve changes between xn and xn1? The Euler method will compute an incorrect value for yn1. This is where the Runge-Kutta methods come into play. The Runge-Kutta methods perform a successive approximation of yn1 by evaluating the f(x,y) function at different locations between xn and xn1. The final computation of yn1 is a linear combination of the successive approximations. For example, the secondorder Runge-Kutta method evaluates f(x,y) at two locations, shown in Eq. (20.7). ¢y1 = ¢xf1xn,yn2 yn + 1 = yn + ¢xfa xn +

1 1 ¢x,yn + ¢y1 b 2 2

(20.7)

The first step of the second-order Runge-Kutta algorithm is the Euler method. A value for y is computed by evaluating f(x,y) at xn. The second step calculates the value of yn1 by evaluating f(x,y) midway between xn and xn1 using a y value halfway between yn and yn  y1. The result is a second-order accurate approximation to yn1. The two-step Runge-Kutta scheme is more accurate than Euler’s method because it does a better job of handling potential changes in slope of f(x,y) between xn and xn1. There are numerous Runge-Kutta schemes of various orders of accuracy. The most commonly used scheme and the one we will implement in this chapter is the fourth-order Runge-Kutta algorithm. As the name implies, it is fourthorder accurate in x. The algorithm consists of five steps, four successive approximations of y and a fifth step that computes yn1 based on a linear combination of the successive approximations. The fourth-order Runge-Kutta solution process is; 1. Find y1 using Euler’s method. 2. Compute y2 by evaluating f(x,y) at a xn +

1 1 ¢x,yn + ¢y1 b . 2 2

3. Calculate y3 by evaluating f(x,y) at axn +

1 1 ¢x,yn + ¢y2 b . 2 2

5568ch20.qxd_jd

278

3/24/03

9:09 AM

Page 278

Chapter 20 Solving Differential Equations

4. Evaluate y4 by evaluating f(x,y) at (xn  x,yn  y3). 5. Compute yn1 using a linear combination of y1 through y4.

The mathematical equations for the five steps are shown in Eq. 20.8. ¢y1 = ¢xf1xn,yn2 ¢y2 = ¢xfaxn +

1 1 ¢x,yn + ¢y1 b 2 2

¢y3 = ¢xfaxn +

1 1 ¢x,yn + ¢y2 b 2 2

(20.8)

¢y4 = ¢xf1xn + ¢x,yn + ¢y32 yn + 1 = yn +

¢y3 ¢y1 ¢y2 ¢y4 + + + 6 3 3 6

Now that we have gone over the derivation of the fourth-order RungeKutta method, let’s write a method to implement it. The method will be named rungeKutta4(). As Runge-Kutta solvers are used by a wide variety of applications, we will define the rungeKutta4() method to be public and static, so it can be universally accessed. We will define rungeKutta4() in a class named ODESolver and place the ODESolver class in the TechJava.MathLib package. The rungeKutta4() method takes three arguments. The first argument is an ODE object (or an ODE subclass object). If you recall, the ODE class will define the number of coupled first-order equations that characterize the ODE and will provide arrays to store the dependent and independent variables. The ODE class also defines the getFunction() method that returns the f(x,y) function for a given x and y. The other two input arguments are the range over which the integration will take place and the increment to the independent variable. This increment will be held constant throughout the entire integration. The number of steps that will be performed is not a user-specified value but is computed based on the range and dx arguments. The integration will stop if the step number reaches the MAX_STEPS parameter defined in the ODE class. The integration follows the steps shown in Eq. (20.8). When the integration is complete, the x[] and y[][] fields of the ODE object will contain the values of the integrated independent and dependent variables. The return value of the rungeKutta4() method is the number of steps computed. The rungeKutta4() code is shown next.

5568ch20.qxd_jd

3/24/03

9:09 AM

Page 279

Runge-Kutta Schemes

package TechJava.MathLib; public class ODESolver { public static int rungeKutta4(ODE ode, double range, double dx) { // //

Define some convenience variables to make the code more readable

int numEqns = ode.getNumEqns(); double x[] = ode.getX(); double y[][] = ode.getY(); //

Define some local variables and arrays

int i,j,k; double scale[] = {1.0, 0.5, 0.5, 1.0}; double dy[][] = new double[4][numEqns]; double ytmp[] = new double[numEqns]; // //

Integrate the ODE over the desired range. Stop if you are going to overflow the matrices

i=1; while( x[i-1] < range && i < ODE.MAX_STEPS-1) { // //

Increment independent variable. Make sure it doesn't exceed the range. x[i] = x[i-1] + dx; if (x[i] > range) { x[i] = range; dx = x[i] - x[i-1]; }

//

First Runge-Kutta step ode.getFunction(x[i-1],dy[0],y[i-1]);

//

Runge-Kutta steps 2-4 for(k=1; k