First let's use the "shooting method". With this method we solve the problem as an initial value problem, and use root finding to match the value on the other side. For example, with

, with and .

we start by solving , with , and , with , such that the first IVP produces and the second IVP produces . Now, using bisection we can adjust until we match .

Let's try it with

with and .

x0 = 1;

v1 = 0;

v2 = 1;

%%% Test v1 and v2

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v1]);

Y(end, 1)

ans = 1

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0, v2]);

Y(end,1)

ans = 4.0000

That's no good because they are both less than 5. Lets not change v2, but lets try to find v1 that will give us a value greater than 5.

x0 = 1;

v1 = 2;

v2 = 1;

xT = 5;

%%% Test v1 and v2

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v1]);

Y(end, 1)

ans = 7.0000

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0, v2]);

Y(end,1)

ans = 4.0000

This works because one is bigger than 5 and the other is less than 5. Now let's write our bisection scheme.

x0 = 1;

v1 = 2;

v2 = 1;

xT=5;

v_mid = (v1 + v2)/2;

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v1]);

t_a = T; % for plotting

x_a = Y(:, 1);

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v2]);

t_b = T; % for plotting

x_b = Y(:, 1);

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v_mid]);

t_mid = T; % for plotting

x_mid = Y(:, 1);

plot(t_a, x_a, t_mid, x_mid, '--', t_b, x_b, 'linewidth', 4)

for i = 1:100

if x_mid(end) == xT

break

else if sign(x_mid(end) - xT) == sign(x_a(end)-xT)

v1 = v_mid;

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v1]);

t_a = T; % for plotting

x_a = Y(:, 1);

else

v2 = v_mid;

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v2]);

t_b = T; % for plotting

x_b = Y(:, 1);

end

end

v_mid = (v1 + v2)/2;

[T, Y] = ode45(@(t, x) myODE(t,x), [0 3], [x0 v_mid]);

t_mid = T; % for plotting

x_mid = Y(:, 1);

plot(t_a, x_a, t_mid, x_mid, '--', t_b, x_b, 'linewidth', 4)

pause(0.1)

end

This very quickly converges to the true solution. To see how the algorithm is working it may be useful to test it out with v1 and v2 that gives us a much larger gap.

We will now start to develop a method for solving boundary value problems numerically. In particular, we will start by solving the problem

, with and .

We already know that the solution is

,

and we would love to find a numerical method that would give us this formula for . Unfortunately, this is generally a hopeless task. (For more complciated equations, it is usually not even possible on paper, and it is certainly not something we can expect MATLAB to do.) Instead, we will try to approximate at many different times. In particular, we will use the same setup that we had for initial value problems and choose a set of evenly spaced t values: , , , ..., , . In general, it only really makes sense to solve a boundary value problem between the two times used for boundary conditions, so we will choose and . Since these times are evenly spaced, we can write each time as , where .

Our goal is to find the values of x at each of these t-values. That is, we would like to find , , , ..., , . Of course, we can't hope to find those values exactly with a numerical method, so we will only obtain approximations , , , ..., , . We therefore need to find an equation (or several equations) for these 's. Actually, we already know two of these points. is just our first boundary condition and is just our second boundary condition. We therefore only have to solve for through . We call these x-values the interior values (as opposed to the boundary values) and the corresponding times through are called the interior points (or interior times). Notice that there are total points, but only interior points.

The key to our approach is to approximate the second derivative in our differential equation with a finite difference scheme. In particular, we already saw in last week's activity that we could approximate the second derivative of a function f with the second order central difference scheme

Translating this into our new notation, we have

This is true for any time t. In particular, we could use the first interior time , which gives us the equation

We don't know what or are, but we do know that this expression is supposed to be approximately equal to zero (because for any time). We can just use this as the definition of our approximations:

Similarly, if we plug in , we get the equation

Again, we know that this is supposed to be approximately equal to zero, so we can take this as a definition for our approximations:

We can repeat this process for any interior time and we get the formula

We therefore have the following system of equations:

,

,

...

,

.

Notice that there are only equations here, because we only plugged in interior times, not the boundary times. (If we used boundary times or in these equations, then we would need x-values from outside of our interval .)

Also note that and are actually already known - they are our boundary values - and so we should really move them over to the right side of our equations. We are then left with

,

...

,

This is a linear system of equations for the interior values , ..., , and so we can rewrite it as a matrix equation . We get

We know several methods to solve a system like this We could use Gaussian elimination, decomposition, or even an iterative method like Jacobi or Gauss-Seidel. We will just use the backslash operator, but this is another example of why linear systems are so important in numerics.

Notice that this method is not a time-stepper, unlike the methods we developed for initial value problems. That is, we can't solve for and then and then , etc. Instead, we have to solve the whole system at once.

As an example, let's implement this in MATLAB to solve

with and .

We will choose .

x0 = 1;

xT = 5;

T = 3;

dt = 0.1;

t = 0:dt:T;

n = length(t);

Notice that n is the total number of points, including the boundary points. There are only interior points, so we need to make an system of equations. We can make a matrix of the right form with the code

v = -2*ones(n-2, 1);

u = ones(n-3, 1);

A = (1/dt^2)*(diag(v) + diag(u, 1) + diag(u, -1));

You already constructed this matrix in a homework problem, so I won't spend any extra time explaining the code above.

We also have to make the vector for the right hand side. This vector is mostly zeros, but the first and last entry depend on the boundary conditions. We can implement this in MATLAB with

b = zeros(n-2, 1);

b(1) = -x0/dt^2;

b(end) = -xT/dt^2;

Finally, we can solve for all of the interior x-values with the command

x_int = A\b;

Keep in mind that this vector does not include the boundary values. Usually we want to know all of the x-values, not just the interior values, but this is easy to fix. The first and last x-value are just the two boudnary conditions, so we can write

x = [x0; x_int; xT];

We would like to plot our solution and compare it to the formula we found by hand. There is a small technical difficulty here. The problem is that the code t = 0:dt:T produces a row vector, but the backslash command produces a column vector. It is unwise to compare vectors of different shapes, so we should make both t and x row vectors or make both t and x column vectors. For no particular reason, let's make t a column vector:

t = t';

Now we can plot our approximation alongside the true solution.

x_true = @(t)((xT - x0)/T * t + x0);

plot(t, x_true(t), 'k', t, x, 'r')

We can't even see both graphs, which is a good sign that we found a good approximation. To confirm, we could check the maximum error

err = max(abs(x - x_true(t)))

err = 1.4655e-14

This is quite small, so we did indeed find a good approximation. In general, we can't expect the error to be quite this small, but it turns out that this differential equation was so simple that our difference scheme was actually exact. Because of that, the only error actually arose from the backslash command and not our BVP solver.

Remember that not all boundary value problems have a solution. In particular, consider

with and

does not have a solution. What happens when we try to solve this problem using our numerical method?

x0 = 0;

xT = 5;

T = pi;

dt = 0.1;

t = (0:dt:T)';

n = length(t);

v = -2*ones(n-2, 1);

u = ones(n-3, 1);

A = (1/dt^2)*(diag(v) + diag(u, 1) + diag(u, -1));

I = eye(n-2);

b = zeros(n-2, 1);

b(1) = -x0/dt^2;

b(end) = -xT/dt^2;

x_int = (A + I)\b;

x = [x0; x_int; xT];

It no longer makes sense to calculate the error, but we can still plot our approximate solution.

plot(t, x, 'r')

We got an answer that looks roughly the same as last time, but remember that there isn't supposed to be a solution. What is going on here? To see the problem, let's try using a smaller . If there was a unique solution, we would expect to get closer to it as we made smaller.

dt = 0.01;

t = (0:dt:T)';

n = length(t);

v = -2*ones(n-2, 1);

u = ones(n-3, 1);

A = (1/dt^2)*(diag(v) + diag(u, 1) + diag(u, -1));

I = eye(n-2);

b = zeros(n-2, 1);

b(1) = -x0/dt^2;

b(end) = -xT/dt^2;

x_int = (A + I)\b;

x = [x0; x_int; xT];

plot(t, x, 'r')

Notice the scale on the y-axis. Our approximations are not getting closer to anything. Instead, they are growing infinitely large as we shrink . This is typical behavior when the BVP doesn't have a solution. The numerical method will always give you an answer, but those answers won't converge as you shrink .

Similarly,

with and has infinitely many solutions. If we try to solve this with our numerical method, we get

x0 = 0;

xT = 0;

T = pi;

dt = 0.1;

t = (0:dt:T)';

n = length(t);

v = -2*ones(n-2, 1);

u = ones(n-3, 1);

A = (1/dt^2)*(diag(v) + diag(u, 1) + diag(u, -1));

I = eye(n-2);

b = zeros(n-2, 1);

b(1) = -x0/dt^2;

b(end) = -xT/dt^2;

x_int = (A + I)\b;

x = [x0; x_int; xT];

plot(t, x, 'r')

In this case, we get , which really is a solution to our differential equation. It turns out that we will still get the same answer no matter how we change . This behavior is a little less common. In general, if your BVP has many solutions, then the numerical solutions might converge to one of those solutions or they might not converge at all.

function dx = myODE(t,x)

dx1 = x(2);

dx2 = 0;

dx = [dx1; dx2];

end