## Topic 14.7: Systems of Initial-Value Problems

Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

We have seen numerous techniques for working with an IVP which consists of a single 1st-order ODE and one initial condition. In this topic, we will look at how we can use the same techniques for solving a system of n 1st-order ODEs with n initial values.

# Background

Useful background for this topic includes:

# References

• Bradie, Section 7.8, Systems of Equations and Higher-Order Equations, p.623.
• Mathews, Section 9.7, Systems of Differential Equations, p.518.

# Theory

Consider the following system of two ODEs:

x(1)(t) = −2 x(t) − 3 y(t) + 3
y(1)(t) = −4 x(t) − 6 y(t) − 2

where we have two initial conditions:

x(0) = −1
y(0) = 5

If we define the vector-valued function u(t) as

and then define

we may rewrite the system of ODEs as:

where f is the vector-valued function

Notice that this form is almost identical to the case with a single ODE:

y(1)(t) = f(t, y(t))
y(t0) = y0

except that all of the functions u(t) and f(t, u) are now vector-valued instead of scalar-valued.

Now, whenever we did any operations in any of the methods, Euler, Huen, Runge-Kutta, or RKF45, the only operations were to:

• Evaluate f(t, y) at some point (t, y) or perhaps (t + h, y + hK),
• Add a scalar multiple of this onto the previous value yk.

These are vector as well as scalar operations, and therefore, we can do the same thing with these functions.

# Example

Letting h = 0.1 and using Euler's method, then

Iterating a number of times, we get the sequence u0 through u6 of:

If we plot the first 21 points u0 through u20, we get the points shown in Figure 1.

Figure 1. The first twenty approximations to the given IVP.

We can solve this IVP exactly to get that:

x(t) = 13/8 e−8 t + 3 t − 21/8
y(t) = 13/4 e−8 t − 2 t + 7/4

If we plot this curve for t between 0 and 2, we get the plot in Figure 2.

Figure 2. The first twenty approximations to the solution of the given IVP.

If x(t) and y(t) represent positions within a plane, then the solution u(t) = (x(t), y(t))T represents the movement of the particle in the plane as time passes. To emphasize this, Figure 3 is an animation showing the path followed by the particle.

Figure 3. An animation of the solution to the IVP.

For this, we used Euler's method, however, it would not be any different to use the 4th-order Runge-Kutta method (or even RKF45). If we use a step size of h = 0, the calculation of u1 would include:

Calculating the weighted average, we get:

Thus, the next point is

This first approximation is shown in Figure 4.

Figure 4. The first approximation using 4th-order Runge-Kutta.

If we plot the first 10 approximations, we get the points in Figure 5.

Figure 5. The first ten approximations using 4th-order Runge-Kutta.

Figure 6 compares the approximations with the actual solution. The approximations are in black while the actual values are in red. The initial significant difference between the approximations and the actual solutions and the closeness of subsequent approximations would suggest that this IVP is well suited for RKF45.

Figure 6. The approximations plotted with the actual solution.

# Problem

Given the vector u(t) = (u1(t), ..., un(t))T, where, for each function uk(t), we have a first-order ODE deined by:

u(1)k(t) = fk(t, u(t))

then we may define

u(1)(t) = f(t, u(t))

where

Given this, together with the initial value u(t0) = u0 = (u0,1, ..., u0,n)T, we would like to approximate the solution to this IVP.

# Assumptions

The function f(t, u) is continous.

# Tools

We will use Euler, Heun, and 4th-order Runge Kutta, though the extension to RKF45 is straight-forward.

# Process

Given a step size h, then tk = t0 + kh. We now examine how each method may be implemented.

### Euler

Given uk, let:

uk + 1 = uk + h f(tk, uk)

### Heun

Given uk, let:

K0 = f(tk, uk)
K1 = f(tk + h, uk + hK0)

Then we define:

### 4th-order Runge Kutta

Given uk, let:

K0 = f(tk, uk)
K1 = f(tk + ½h, uk + ½hK0)
K2 = f(tk + ½h, uk + ½hK1)
K3 = f(tk + h, uk + hK2)

Then we define:

# Example 1

Almost all students have seen the butterfly effect, shown here on Wikipedia.

This is the result of a system of three non-linear differential equations collectively known as the Lorenz equations:

where

These variables represent certain properties of the earth's atmosphere. If we start with an initial condition of the atmosphere, then the differential equation dictates how the properties change over time.

For demonstration purposes, we will use σ = 10, ρ = 28, and β = 2.7

If we start with the initial point u0 = (1, 1, 1)T and use the function

together with a step size of h = 0.01, then applying Euler's method, we can calculate u1:

If we repeat this 1000 times, we can plot the result, as shown in Figure 1.

Figure 1. The first 1000 approximations for the solution to the Lorenz equations (click to enlarge).

This indicates that each of the three values varies over time, however, the first two variables do not move outside the range [-20, 20] while the third variable appears to be bounded by [0, 50].

It is more common to interpolate the points, as is shown in Figure 2.

Figure 2. The first 10000 approximations for the solution to the Lorenz equations (click to enlarge).

That the solution to this differential equation follows such a peculiar path is part of the study of chaos theory.

To be completed.

# Error

The error analysis for a system of IVPs is similar to that of an IVP with just one ODE.

# Question 1

Consider the system of IVPs defined by

x(1)(t) = -0.6 x(t) + 1.8 y(t)
y(1)(t) = -2.6 x(t) − 0.3 y(t)

together with

x(0) = 1 y(0) = 1

Find the first five iterations using h = 0.1.

```>> M = [-0.6 1.8; -2.6 -0.3];
>> x = zeros( 2, 1000 );
>> x(:, 1) = [1 1]';
>> for i=2:1000
x(:, i) = x(:, i - 1) + 0.1 * M*x(:, i - 1);
end;
>> plot( x(1, :), x(2, :) )
```

however, the first two iterations are u1 = (1.12, 0.71)T and u2 = (1.1806, 0.3975)T.

# Matlab

If the IVP is linear, that is, the function f may be written as

f(t, u) = Mu + b

for a matrix M and vector b, then we can very easily use Matlab to solve such systems.

For example, given the IVP

x(1)(t) = -0.5 x(t) + 1.9 y(t) + 0.4
y(1)(t) = -2.6 x(t) − 0.7 y(t) + 0.8
x(0) = 1
y(0) = 1

then we may assign M and b:

```>> M = [-0.5 1.9; -2.6 0.7];
>> b = [0.4 0.8]';
```

We will store the iterations in the 2 × N matrix and assign the initial values to be the first column:

```>> N = 10000;
>> x = zeros( 2, N );
>> x(:, 1) = [1 1]';
```

Next, we iterate Euler's method by assigning to the next column:

```>> h = 0.01;
>> for i=2:N
x(:, i) = x(:, i - 1) + h*(M*x(:, i - 1) + b);
end;
```

To plot this data, we pass plot the first row (the x values) and the second row (the y values):

```>> plot( x(1, :), x(2, :) )
```

The result of the first 10000 iterations is shown in Figure 1.

Figure 1. The first 10000 approximations of the solution to the system of IVPs (click to enlarge).

You may note that the solution appears to converge to a point. That point is the value of u such that Mu + b = 0, or Mu = −b. Thus, to find this point, we solve:

```>> M \ (-b)
ans =

0.340264650283554
-0.120982986767486
```

# Heun's Method

To implement Heun's method, we need only one change:

```>> for i=2:N
K0 = M*x(:, i - 1) + b;
K1 = M*(x(:, i - 1) + h*K0) + b;
x(:, i) = x(:, i - 1) + h*(K0 + K1)/2;
end;
```

# 4th-order Runge Kutta Method

To implement the 4th-order Runge Kutta method, we need only one change:

```>> for i=2:N
K0 = M*x(:, i - 1) + b;
K1 = M*(x(:, i - 1) + h/2*K0) + b;
K2 = M*(x(:, i - 1) + h/2*K1) + b;
K3 = M*(x(:, i - 1) + h*K2) + b;
x(:, i) = x(:, i - 1) + h*(K0 + 2*K1 + 2*K2 + K3)/6;
end;
```

# Maple

Maple works with an IVP consisting of a systems of ODEs as easily as it does with a single ODE.

```dsolve( {
D(x)(t) = -2*x(t) − y(t) + 2,
D(y)(t) = −x(t) − 2*y(t) − 1,
x(0) = 1,
y(0) = 1
} );
```

This results in the solution

{x(t) = 5/6 e − 3 t−3/2 et + 5/3, y(t) = 5/6 e−3 t+3/2 et − 4/3}
```dsolve( {
D(x)(t) = -0.5*x(t) + 1.9*y(t) + 0.4,
D(y)(t) = -2.6*x(t) - 0.7*y(t) + 0.8,
x(0) = 1,
y(0) = 1
} ):
evalf(%);

{ x(t) =  0.34026 + exp(-0.6 t) ( 0.98896 sin(2.2204 t) + 0.65973 cos(2.2204 t)),
y(t) = -0.12098 + exp(-0.6 t) (-0.82302 sin(2.2204*t) + 1.12098 cos(2.2204 t)) }
```

Plotting this solution for t = [0, 15], we get the result shown in Figure 1.

Figure 1. The solution to the initial value problem shown above.