Paul Heckbert

2 Nov. 2000, revised 17 Dec. 2000

I illustrate shooting methods, finite difference methods, and the collocation and Galerkin finite element methods to solve a particular ordinary differential equation boundary value problem.

- y''=-g (the differential equation)
- g=9.8 [m/sec^2] (acceleration due to gravity)
- y(0)=0 (launch from ground)
- y(5)=40 (fireworks explode after 5 seconds, we want them 40 m off ground)

This can be solved analytically; exact motion is quadratic in t, final answer is y'(0)=32.5 [m/sec].

h=.5 - h too big

h=.1 - smaller h gives more accurate results.
But note that the y'(0) that secant method solves for, in red,
is still not correct
(not 32.5), because of errors of our IVP solution.
If we used Runge-Kutta or other more accurate IVP solution method, instead
of Euler's method, that would help.

Below is my Matlab code for the shooting method.

function y = rocket(dy0) % return altitude at t=te (y(te)) as a function of initial velocity (y'(0)=dy0) global h te [tv,yv] = euler2(h,0,te,0,dy0); plot(tv,yv,'o-', 'LineWidth',2); % invariant: tv(te/h+1)==te y = yv(te/h+1); % returns y at t=te return; function [tv,yv] = euler2(h,t0,tmax,y0,dy0) % use Euler's method to solve 2nd order ODE y''=-g+a*y^2 % return tvec and yvec sampled at t=(t0:h:tmax) as col. vectors % y(t0)=y0, y'(t0)=dy0 % Paul Heckbert 30 Oct 2000 global a; % coeff of nonlinear acceleration g = 9.8; % accel. of gravity, [m/sec^2] y = y0; % position dy = dy0; % velocity tv = [t0]; yv = [y0]; for t = t0:h:tmax y = y+dy*h; % this and following line are Euler's method dy = dy+(-g+a*y^2)*h; tv = [tv; t+h]; yv = [yv; y]; end return;

function ret = shoot(hh) % shooting method for fireworks problem global a te ye h; a = 0; % coeff of nonlinear acceleration te = 5; % end time [sec] ye = 40; % end height [m] h = hh; clf; hold on; for dy = 20:10:50 y = rocket(dy); text(te+.2,y,sprintf('y\047(0)=%g',dy), 'FontSize',15); end % now use root finding to find correct initial velocity dy = secant(20,30,1e-4); % draw last curve y = rocket(dy); text(te+.2,y,sprintf('y\047(0)=%g',dy), 'Color','r', 'FontSize',15); set(gca, 'FontSize', 16); % for tick marks line([te te],[-40 140], 'Color','k'); line([te te],[ye ye], 'Color','r', 'Marker','o', 'LineWidth',3); xlabel('t', 'FontSize',20); ylabel('y', 'FontSize',20); title(sprintf('Shooting Method on y\047\047=-g'), 'FontSize',20); return; function x = secant(x1,x2,tol) % secant method for one-dimensional root finding global ye; y1 = rocket(x1)-ye; y2 = rocket(x2)-ye; while abs(x2-x1)>tol disp(sprintf('(%g,%g) (%g,%g)', x1, y1, x2, y2)); x3 = x2-y2*(x2-x1)/(y2-y1); y3 = rocket(x3)-ye; x1 = x2; y1 = y2; x2 = x3; y2 = y3; end x = x2; return;

When we do shooting method we find roots of this function.
For this particular, very simple, problem the function is linear,
so root finding is trivial (secant method converges in one iteration).
In general, the root-finding problem is linear if the
differential equation is linear.

Above plot is made like so:

function rocketplot % plot the function which we're finding roots of, in shooting method global a; a = 0; % coeff of nonlinear acceleration hold off; yv = []; dy0v = 20:5:50; for dy0 = dy0v y = rocket(dy0); yv = [yv y]; end set(gca, 'FontSize', 16); % for tick marks plot(dy0v,yv, '-', 'LineWidth',2); xlabel('initial velocity', 'FontSize',20); ylabel('y(5) = altitude at t=5', 'FontSize',20); title(sprintf('Root-Finding Problem for Shooting Method on y\047\047=-g'), 'FontSize',16); line([20 50], [40 40], 'Color','k');

Now let's temporarily look at a different problem, y''=-g+.01*y^2, which is nonlinear. (I know of no reasonable physical basis for this ODE; it's cooked up.)

Trajectories are no longer parabolas.

And makes this function nonlinear, so root finder has to do more work
(takes about 5 iterations).

This ODE was solved with code just like shoot.m, shown earlier, except a=.01 .

Now we go back to original problem, y''=-g.

With n=10 intervals and n+1=11 function samples.
We end up solving an 11x11 tridiagonal system.
Can be done in time linear in n.

function findif(n) % finite differences method to solve fireworks problem % n is number of divisions global g te ye h; g = 9.8; % accel. of gravity, [m/sec^2] te = 5; % end time [sec] ye = 40; % end height [m] h = te/n; A(1,1) = 1; b(1) = 0; % y0 = 0 for i = 2:n A(i, i-1:i+1) = [1 -2 1]; b(i) = -g*h^2; end A(n+1,n+1) = 1; b(n+1) = ye; b = b'; % make b into a column vector y = A\b; A b' y' (y(2)-y(1))/h (4*y(2)-y(3))/(2*h) % best est. of initial velocity (parab thru 1st 3 pts) hold off; set(gca, 'FontSize', 16); % for tick marks plot(0:h:te, y, 'bo-', 'LineWidth',2); xlabel('t', 'FontSize',20); ylabel('y', 'FontSize',20); title('Finite Difference Method', 'FontSize',20);Here is what the above code prints out.

>> findif(10) A = 1 0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0 1 b = Columns 1 through 7 0 -2.4500 -2.4500 -2.4500 -2.4500 -2.4500 -2.4500 Columns 8 through 11 -2.4500 -2.4500 -2.4500 40.0000 y = Columns 1 through 7 0 15.0250 27.6000 37.7250 45.4000 50.6250 53.4000 Columns 8 through 11 53.7250 51.6000 47.0250 40.0000 ans = 30.0500 linear est. of initial velocity ans = 32.5000 quadratic est. of initial velocity

Note that the second estimate of launch velocity above is exact.

Blue: solution u(t); Red: scaled basis functions u[i]*phi[i](t) for i=1..n+2.
With n=4 interior intervals and n+2=6 basis functions.
We end up solving a 6x6 tridiagonal system.
Can be done in time linear in n.

Because the correct solution is quadratic, and our basis functions span the space of all quadratic functions, collocation computes an exact solution to this problem, even when n=1.

function colloc(n) % collocation finite element method to solve fireworks problem % n is number of divisions % using quadratic B-spline basis functions global g te ye h; g = 9.8; % accel. of gravity, [m/sec^2] te = 5; % end time [sec] ye = 40; % end height [m] h = te/n; A(1, 1:2) = [1/2 1/2]; b(1) = 0; % y=0 at t=0 for i = 2:n+1 A(i, i-1:i+1) = [1 -2 1]; b(i) = -g*h^2; end A(n+2, n+1:n+2) = [1/2 1/2]; b(n+2) = ye; % y=ye at t=te b = b'; % make b into a column vector u = A\b; A b' u' (u(2)-u(1))/h % initial velocity % plot(-h/2:h:te+h/2, u, 'g'); % plot the coefficients tv = 0:h/8:te; % supersample the parabolas yv = []; for t = tv i = floor(t/h)+2; if i==n+2 i = n+1; end x = t/h-i+1.5; y = u(i-1)*.5*(x-.5)^2 + u(i)*(.75-x^2) + u(i+1)*.5*(x+.5)^2; yv = [yv y]; end % compute joints between elements tk = (0:h:te)'; yk = .5*(u(1:n+1)+u(2:n+2)); clf; hold on; set(gca, 'FontSize', 16); % for tick marks plot(tv,yv,'b', tk,yk,'bo', 'LineWidth',2); xlabel('t', 'FontSize',20); ylabel('y', 'FontSize',20); title('Collocation Finite Element Method', 'FontSize',20); % draw scaled basis functions dh = h/15; eps = 1e-5; for i = 1:n+2 tv1 = (i-3)*h:dh:(i-2)*h-eps; yv1 = u(i)*.5*(tv1/h-i+3).^2; tv2 = (i-2)*h:dh:(i-1)*h-eps; yv2 = u(i)*(.75-(tv2/h-i+1.5).^2); tv3 = (i-1)*h:dh:i*h+eps; yv3 = u(i)*.5*(tv3/h-i).^2; plot([tv1 tv2 tv3], [yv1 yv2 yv3], 'r', 'LineWidth',2); end hold off;Here is what the above code prints out:

>> colloc(4) A = 0.5000 0.5000 0 0 0 0 1.0000 -2.0000 1.0000 0 0 0 0 1.0000 -2.0000 1.0000 0 0 0 0 1.0000 -2.0000 1.0000 0 0 0 0 1.0000 -2.0000 1.0000 0 0 0 0 0.5000 0.5000 b = 0 -15.3125 -15.3125 -15.3125 -15.3125 40.0000 u = -20.3125 20.3125 45.6250 55.6250 50.3125 29.6875 ans = 32.5000

Blue: solution u(t); Red: scaled basis functions u[i]*phi[i](t) for i=1..n+2.
With n=4 interior intervals and n+2=6 basis functions.
We end up solving a 6x6 pentadiagonal system.
Can be done in time linear in n.

Note that solution is exactly the same as with collocation method, for this particular problem, and is an exact solution, but that's not true in general. We could also solve this problem exactly with n=1 using these basis functions.

function galerkin(n) % Galerkin finite element method to solve fireworks problem % n is number of divisions % using quadratic B-spline basis functions global g te ye h; g = 9.8; % accel. of gravity, [m/sec^2] te = 5; % end time [sec] ye = 40; % end height [m] h = te/n; A(1, 1:2) = [.5 .5]; b(1) = 0; % y=0 at t=0 A(2, 1:4) = [4 -7 2 1]; % Galerkin condition on first interior intvl b(2) = -5*g*h^2; for i = 3:n A(i, i-2:i+2) = [1 2 -6 2 1]; b(i) = -6*g*h^2; end A(n+1, n-1:n+2) = [1 2 -7 4]; % last interior interval b(n+1) = -5*g*h^2; A(n+2, n+1:n+2) = [.5 .5]; b(n+2) = ye; % y=ye at t=te b = b'; % make b into a column vector u = A\b; A b' u' (u(2)-u(1))/h % initial velocity % plot(-h/2:h:te+h/2, u, 'g'); % plot the coefficients tv = 0:h/8:te; % supersample the parabolas yv = []; for t = tv i = floor(t/h)+2; if i==n+2 i = n+1; end x = t/h-i+1.5; y = u(i-1)*.5*(x-.5)^2 + u(i)*(.75-x^2) + u(i+1)*.5*(x+.5)^2; yv = [yv y]; end % compute joints between elements tk = (0:h:te)'; yk = .5*(u(1:n+1)+u(2:n+2)); hold off; set(gca, 'FontSize', 16); % for tick marks plot(tv,yv,'b', tk,yk,'bo', 'LineWidth',2); xlabel('t', 'FontSize',20); ylabel('y', 'FontSize',20); axis([0 5 0 60]); title('Galerkin Finite Element Method', 'FontSize',20);Here is what the above code prints out.

>> galerkin(4) A = 0.5000 0.5000 0 0 0 0 4.0000 -7.0000 2.0000 1.0000 0 0 1.0000 2.0000 -6.0000 2.0000 1.0000 0 0 1.0000 2.0000 -6.0000 2.0000 1.0000 0 0 1.0000 2.0000 -7.0000 4.0000 0 0 0 0 0.5000 0.5000 b = 0 -76.5625 -91.8750 -91.8750 -76.5625 40.0000 u = -20.3125 20.3125 45.6250 55.6250 50.3125 29.6875 ans = 32.5000Sometimes you can get the same accuracy with a smaller n, using Galerkin instead of collocation. So although it's more algebra and calculus to set up the Galerkin system, and more computation to solve it, it is often worth the trouble to use Galerkin.

We've discussed three methods: shooting, finite difference, and finite element. All of these methods transform boundary value problems into algebraic equation problems (a.k.a. root-finding). When the differential equation is linear, the system of equations is linear, for any of these methods. When the differential equation is nonlinear, the system of equations is, in general, nonlinear. Broyden's method, or other nonlinear system solver, could be used in the latter case.

Not discussed here: superposition method. Very useful for homogeneous linear ODE's.

15-859B, Introduction to Scientific Computing