This month, we will look at the format required to solve a differential equation or a system of differential equations using one of the command-line differential equation solvers such as rkfixed, Rkadapt, Radau, Stiffb, Stiffr or Bulstoer.

For a numerical routine to solve a differential equation (DE), the differential equation must be passed as an argument to the solver routine. A standard form for all DEs allows you to do this.

Here's the basic idea: given a differential equation, get rid of any second, third, fourth, and so on derivatives that appear, leaving only first derivatives.

To see what this means, look at a fairly simple differential equation.

This DE contains a second derivative. How do you write a second derivative as a first derivative? A second derivative is a first derivative of a first derivative.

Define two functions y₀(x) and y₁(x) as

Then your differential equation can be written as

Now this DE contains two functions instead of one, but there is a strong relationship between these two functions, namely,

Your original DE is now a system of two DEs,

Further, adopt the convention of writing these equations with the derivatives alone on the left-hand side.

This is the first step in the standardization process. Below are additional examples of converting a DE into a system involving only first order derivatives.

DE 1

The differential equation

has a system representation of

where

is the same as

Note that in this example the naming convention has been such that

for each n; that is, the subscript indicates the order of the derivative of the original unknown function. This is a useful convention, which you should always use.

DE 2

The differential equation

corresponds to a system (where U₀(r) = U(r))

Note that to get this last equation, you had to solve for

by collecting terms on one side and using sin to get rid of the asin.

The next step needed for the numerical solvers is to capture the information in this system in a single function. Return to the first system derived above.

There are two functions of x,

and

Suppose at a particular value of x,

you know the values of these two functions

Put these values in a vector y.

Given the values of the function at a value x, what are the values of

and

at

_?

The DE system tells you that

At each x, the derivative of y₀ is equal to the value of y₁, and the derivative of y₁ is nicely expressed in terms of the values of y₀ and y₁. Thus, the derivatives of y₀ and y₁, respectively, at

are

There are general formulas as well. The function

specifies for each vector of function values at the point x, the corresponding vector of derivative values.

This function DY is the essential component to be used in solving a differential equation numerically. When constructing this function, be sure to use the array subscript to indicate the order of the derivative.

Construct the corresponding function for the example systems from above.

DE 1

Define

DE 2

Define

Of course, this whole process is described in great detail. After you've constructed the function D a few times, you'll be able to do it quite easily for a given ODE.

For instance, look at the following DE and the corresponding function. See how much easier it is to construct the function D:

Systems of Differential Equations

It may be that you want to solve a system of equations rather than a single differential equation. Consider the system

You can define

and then the system looks like

By getting each derivative alone on one side, you have

which is the standard form for a system. You can now write the associated function DF:

It isn't that much different from what you did earlier, but in this case you have to remember what each element F_i stands for. In this case both F₀ and F₁ are function values, not derivatives.

For a more complicated example, consider a second-order system:

Define

Then, the system of two equations becomes the system

and the corresponding function is

In this case, the functions are defined so that the even indices correspond to the unknown function u and its derivatives, while the odd indices correspond to v.

Now that you've set up the derivative vector, solve the second-order system using Rkadapt.

Rkadapt takes the following arguments:

Rkadapt(ic, start, end, npts, D)

• ic - vector of initial conditions. Scalar for a first-order DE.
• start, end - endpoints of the solution interval
• npts - number of approximating points past the initial condition. The result is a matrix containing npts + 1 rows
• D - the derivative vector

Initial condition vector

Vector-valued function containing the derivatives of the unknown functions

Rkadapt returns a vector of data where the first column contains the values of the independent variable. The remaining columns contain the function and derivatives evaluated at the corresponding values of the independent variable. For a system, the columns are returned in the same order in which the initial conditions are entered.

Compare Results to Odesolve

Was this article helpful? Let us know.