Boundary Value Problems

15-859B, Introduction to Scientific Computing
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.


The Problem

Solve for y(t), altitude of rocket, given In particular, what should launch velocity y'(0) be?

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


Shooting Method

Matlab code for this 2nd order ODE using Euler's method:


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.

Listing of rocket.m

    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;

Listing of shoot.m

    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:

Listing of rocketplot.m

    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.


Finite Difference Method

We let y[1+i] = y(i*h) and discretize the problem.


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.

Listing of findif.m

    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.

Listing of findif.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.


Collocation Finite Element Method

Because we need 2nd derivatives, we should use basis functions that are at least quadratic in t. I chose quadratic B-spline basis functions. (Note that many other choices are possible here.) Approximate solution is a linear combination of these u(t) = sum{i=1..n+2}{u[i]*phi[i](t)}. We constrain the residual r(t)=u''(t)+g to be zero at the center of each interior interval.


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.

Listing of colloc.m

    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:

Listing of colloc.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

Galerkin Finite Element Method

In the Galerkin method, the residual is constrained to be orthogonal to each of the basis functions phi[i](t), for the interior intervals i=2..n+1.


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.

Listing of galerkin.m

    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.

Listing of galerkin.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.5000
Sometimes 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.

Discussion

If you wanted to compare the above methods, you could plot error as a function of flops for each, and see which wins.

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