Week 9 Lecture 3: Boundary Value Problems
Rahman notes:
Shooting Method
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 
.
, 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.
Finite Differences
I borrowed the old notes for this, but changed a few things up in the video.
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 
.
. 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.
, , , ..., , . 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
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
, 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:
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
, 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
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 
.)
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
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
, ..., , 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.
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.
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
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
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.
Problems With No Solution
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.
. 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 .
. 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 . Problems With Many Solutions
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.
, 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