Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

In this topic, we solve one-dimensional boundary-value problems.

Theory

Given the linear 2nd-order ordinary differential equation

c₂(t) y⁽²⁾(t) + c₁(t) y⁽¹⁾(t) + c₀(t) y(t) = g(t)

where the c_i(t) and g(t) are known functions, a boundary-value problem is a problem where the solution y(t) is subject to two constraints, each given at a different value of t:

y(a) = y_a
y(b) = y_b

Thus, at each interior point on [a, b], the solution y(t) must satisfy the differential equation. We can break the interval [a, b] into n sub-intervals by defining h = (b − a)/n and letting t_i = a + ih. Therefore, we see that:

t₀ = a
t_n = b

We would like to find approximations of y(t_i), so let us represent the approximation of this by y_i ≈ y(t_i). From the initial conditions, we already have two such values:

y₀ = y_a
y_n = y_b

Recall from Topic 13 (Differentiation) that we may approximate the first and second derivatives by the centred divided-difference formulae:

If we substitute these two approximations into our ODE, we get its a discrete form, that is, a difference equation:

If we multiply through by 2h² and collect with respect to the three values y_{i + 1}, y_i, and y_{i − 1}, we get the form:

If we evaluate this for each value of i = 1, 2, ..., n − 1, we get n − 1 equations:

Notice that in the first and last equations, y₀ and y_n are known. Therefore, we have n − 1 unknowns, namely, y₁, ..., y_{n − 1}.

We may solve this system of linear equations using any of the tools we have available from linear algebra.

Special Case: Constant Coefficients

In the special case where the ODE has constant coefficients, this greatly simplifies the definition of the system of linear equations. If the matrix is defined as

c₂ y⁽²⁾(t) + c₁ y⁽¹⁾(t) + c₀ y(2)(t) = g

Writing out the equations, we note that all of the coefficients are identical:

Therefore, if we define the four variables l, d, u, and v (for lower, diagonal, upper, and vector, respectively):

l = 2c₂ − hc₁
d = 2h²c₀ − 4c₂
u = 2c₂ + hc₁
v = 2h²g

We can rewrite the equations as:

       dy₁ + uy₂ = v − ly₀
ly₁ + dy₂ + uy₃ = v
ly₂ + dy₃ + uy₄ = v
              .
              .
              .
ly_{n − 2} + dy_{n − 1} = v − uy_n

We can write this as the system of linear equations:

which we can easily solve for y.

Once we have the vector y, we have approximations for each of the interior points y₁, ..., y_{n − 1}. If we needed a point between these, we could use interpolation.

HOWTO

Problem

Approximate the solution of the 2nd-order linear ordinary differential equation

c₂(t) y⁽²⁾(t) + c₁(t) y⁽¹⁾(t) + c₀(t) y(2)(t) = g(t)

subject to the boundary conditions

y(a) = y_a
y(b) = y_b

Assumptions

The functions c_i(t) are appropriately continuous.

Tools

We will use the centred divided-difference formula for the derivative and linear algebra.

Process

Divide the interval [a, b] into n sub-intervals by defining h = (b − a)/n, setting t_i = a + ih for i = 0, 1, 2, ..., n. Let y_i represent the approximation of y(t_i), and therefore

y₀ = y_a
y_n = y_b

Rewrite the differential equation

We can simplify this by multiplying by 2h² and collecting on the y's:

Evaluate this equation at each of the points i = 1, 2, ..., n − 1.

This defines a system of n − 1 linear equations and n − 1 unknowns. Thus, this can be written in the form My = g.

Special Case

In the special case where the coefficients are constant, that is, the differential equation is of the form:

c₂ y⁽²⁾(t) + c₁ y⁽¹⁾(t) + c₀ y(2)(t) = g

and therefore the ODE simplifies to

Now define:

l = 2c₂ − hc₁
d = 2h²c₀ − 4c₂
u = 2c₂ + hc₁
v = 2h²g

and solve the system of linear equations:

Examples

Example 1

Solve the boundary-value problem

y⁽²⁾(t) + 3 y⁽¹⁾(t) + 8 y(t) = 0

subject to

y(0) = 1
y(1) = 2

at 9 interior points.

Using n = 10 and therefore h = 0.1, we can find:

l = 2 − 0.3 = 1.7
d = 0.16 − 4 = -3.84
u = 2 + 0.3 = 2.3
v = 0

Thus, we are solving the system

Solving this yields

If we plot these points and the actual solution (y(t) ≈ 6.6199 e^−1.5 t(2.1642 sin(2.3979 t) + 0.1511 cos(2.3979 t))) we get plot shown in Figure 1.

Figure 1. The solution to the BVP for Example 1 together with the approximation.

If we wanted a better approximation, we could use a smaller value of h.

Engineering

To be completed.

Error

1st-Order Boundary-Value Problems

To determine the error for the 1st-order backward divided-difference formula, we need only look at the Taylor series approximation:

Simply rearranging and dividing by h yields the formula:

Thus, the error is O(h).

2nd-Order Boundary-Value Problems

To determine the error for the 2nd-order backward divided-difference formula, we need only look at the two Taylor series approximations:

and subtract the first from 4× the second, and rearrange, we get the that:

Note that the sum of the coefficients in the parentheses is 1, and therefore may be approximated by the average of g⁽³⁾(x) on the interval [x - 2h, x], and therefore the formula is O(h²).

Compare the error of the 2nd-order formula to that using the error of the 2nd-order centred divided-difference formula which has a coefficient -1/6, and thus, the centred divided-difference formula has, approximately, half the absolute error of the backward divided-difference formula. To view this, consider Figure 1. Here we see the points used to approximate the 2nd-order backward (black) and centred (blue) divided-difference formulae to approximate the derivative (magenta) at the fixed point. You will note that the error of the centred approximation is approximately half that of the backward approximation.

Figure 1. Comparison of 2nd-order centred and backward divided-difference approximations of the derivative.

The function shown in Figure 1 is g(x) = exp( x² ) and the point is x = 0.5 . Table 1 shows the approximations and the errors for h = 0.1 and h = 0.01. The actual derivative at that point is 1.284025417.

Table 1. Comparison of errors.

h	Centred	Centred error	Backward	Backward error
0.1	1.2991	0.0151	1.2610	0.023
0.01	1.2842	0.00015	1.2837	0.00029

You will notice that the error of the centred divided-difference formula is approximately half that of the backward divided-difference formula. Additionally, both errors in the second row are approximately 0.1² = 0.01 that of the first row.

Questions

Question 1

Solve the boundary value problem

y⁽²⁾(t) + 2 y(t) = 1 y(0) = 0
y(2) = 1

using h = 0.5.

Answer: (0.5, 0.16667), (1.0, 0.5), (1.5, 0.83333)

Matlab

Maple

Given the constant-coefficient boundary value problem

c₂ y⁽²⁾(t) + c₁ y⁽¹⁾(t) + c₀ y(2)(t) = f y(a) = y_a
y(b) = y_b

we can solve such a system in Maple as follows:

c[2] := 1.0;
c[1] := 4.0;
c[0] := 2.0;
f := 0.0;
n := 10;
a := 1.0;
b := 3.0;
ya := 4.0;
yb := 2.0;
h := (b - a) / n;
low := 2*c[2] - h*c[1];
diag := 2*h^2*c[0] - 4*c[2];
up := 2*c[2] + h*c[1];
vec := 2*h^2*f;
M := Matrix( n - 1, n - 1 );
v := Vector( n - 1 );
M[1, 1] := diag;
M[1, 2] := up;
v[1] := vec - ya*low;

for i from 2 to n - 2 do
    M[i, i - 1] := low;
    M[i, i] := diag;
    M[i, i + 1] := up;
    v[i] := vec;
end do:

M[n - 1, n - 2] := low;
M[n - 1, n - 1] := diag;
v[n - 1] := vec - yb*up;

y := LinearAlgebra[LinearSolve]( M, v );

Topic 15.1: Finite-Difference Method