Two-point BVPs for ODEs

• We consider two methods for solving two-point boundary value problems using the shooting method and finite difference method.

• The example problem is given as a scalar second-order ODE,

  

  with boundary conditions

  

The shooting method

• Please read the lecture note 11.

An example Matlab code: Click to download the sample code bvp_shooting.m rootFind.m

% --------------------------------------------%
% The 2-point BVP for ODEs -- shooting method %
%         Written by Prof. Dongwook Lee       %
%           AMS 213B, Spring 2016             %
% --------------------------------------------%

% --------------------------------------------%
% BVP:
% A 2nd order ODE: u''=6t, 0<t<1, 
% together with BCs: u(0)=0 & u(1) = 1
%
% Convert it to a system of 1st order ODEs:
%
% [y1]   [y2]
% |  | = |  |
% [y2]   [6t],
%
% where y1=u, y2=y1'=u'
%
% BVP ==> Use shooting method to solve the IVP 
% with y1(ta)=alpha, y2(ta)=beta
% --------------------------------------------%

clear all;
format long;

% temporal discretization
tStep = 128;

% two end boundary points
ta=0;
tb=1;

% time step
dt=(tb-ta)/tStep;

% set an array for time t
tArr=linspace(ta,tb,tStep+1);
y1(tStep)=0;
y2(tStep)=0;

% set two given BCs: u(0) = 0 & u(1) = 1
ua=0;
ub=1;

% set initial threshold value to be large
epsilon=100;

% set counter
nCounter=1;

% set the very first initial guess for the slope u'(0) for shooting
y2_ta_init = 3;


% time marching with two exit conditions:
%   exit when either epsilon is smaller than 1.e-7
%   or the number of time marching is more than 100.
while (epsilon > 1.e-7) & (nCounter < 100)
    
    % set clock to zero
    t=0;
    
    % set the left BC
    y1(1)=ua;

    % the inner while loop implements the shooting method with
    % initial guesses
    while (t<=tb);
        if (nCounter == 1)
            y2(1)=y2_ta_init; %initial guess for slope u'(0)
        end
        
        % let's store all attempted initial guesses
        y2_guess_at_a(nCounter)=y2(1);

        % Forward Euler for y1 and y2
        for i=1:tStep;
            y1(i+1)=y1(i)+dt*y2(i);
            y2(i+1)=y2(i)+dt*6*(i-1)*dt;
            %pause
        end
        
        % update time
        t=t+dt;
    end
    
    % calculate the error
    y1_error(nCounter)=y1(end)-ub;
    
    % Call the root finding routine: by calling the root finder,
    % we keep updating the initial guesses y2(1), with which
    % we reiterate the shooting (the inner while loop)
    [y2(1),epsilon]=rootFind(y2(1),y1(end)-ub,y2(end));
    
    %epsilon;
    
    % set counter
    nCounter=nCounter+1;
    
    % plot all iterated solutions as a movie
    figure(1);plot(tArr,y1);pause(0.1);
end

% show the total number of steps it takes for the final solution
nCounter

% plot the final converged solution
figure(1);plot(tArr,y1,'k*-');

% plot the history of initial guesses
figure(2);plot(y2_guess_at_a,'r*-');

% plot the history of y1 errors
figure(3);plot(y1_error,'b*-');

% --------------------------------------------%
% The 2-point BVP for ODEs -- shooting method %
%         Written by Prof. Dongwook Lee       %
%           AMS 213B, Spring 2016             %
% --------------------------------------------%
% Root finding using the exact h' = y_2(tb)


function [y,epsilon] = rootFind(y,h,hprime);

    epsilon = h/hprime;
    y=y-epsilon;
    epsilon=abs(epsilon);

Note

• Depending on what initial guess (i.e., the initial slope of the shooting) you choose the solution of the shooting method may not converge, or if it converges it does so with different convergence rates (i.e., faster convergence with a better initial guess).

The finite difference method

• Please read the lecture note 11.

An example Matlab code: Click to download the sample code bvp_fd.m

% --------------------------------------------%
% The 2-point BVP for ODEs: finite diff method%
%         Written by Prof. Dongwook Lee       %
%           AMS 213B, Spring 2016             %
% --------------------------------------------%

% --------------------------------------------%
% BVP:
% A 2nd order ODE: u''=6t, 0<t<1, 
% together with BCs: u(0)=0 & u(1) = 1
%
% Convert it to a system of 1st order ODEs:
%
% [y1]   [y2]
% |  | = |  |
% [y2]   [6t],
%
% where y1=u, y2=y1'=u'
%
% BVP ==> Use a finite difference method to solve the IVP 
% with y1(ta)=alpha, y2(ta)=beta
% --------------------------------------------%

%BVP u''=6t, 0<t<1, u(0)=0 & u(1) = 1
%This routine implements a finite difference method

clear all;
format long;

% temporal discretization
tStep = 8;

% two end boundary points
ta=0;
tb=1;

% time step
dt=(tb-ta)/tStep;

% set an array for time t
tArr=linspace(ta,tb,tStep+1);

% set two given BCs: u(0) = 0 & u(1) = 1
ua=0;
ub=1;

% set up initial arrays of the matrix A and 
% the vector b for Ax = b
a(1:tStep-1,1:tStep-1)=0;
b(1:tStep-1)=0;

% initialize A
for j=1:tStep-1
    for i=1:tStep-1
        if (abs(i-j)<=1)
            if (i ~= j) 
                a(i,j)=1;
            else
                a(i,j)=-2;
            end
        end
    end
end
% make sure we rescale with dt^2
a=a/dt^2;


% initialize b 
for j=1:tStep-1;
    b(j)=6*j*dt;
end

% make sure we incorporate the BCs
b(1)=b(1)-ua/dt^2;
b(end)=b(end)-ub/dt^2;

% solve Ax = b by direct inversion
x=inv(a)*transpose(b);

% plot the final solution
figure(4);plot(tArr,[ua;x;ub],'ko-')

Note

• The finite difference approach converts the BVP into a linear system .
• It is crucial to include the two boundary conditions in the right hand side vector .
• One needs to solve the linear system using an efficient numerical linear algebra solver to invert . Although the example MATLAB code directly inverts the matrix, in practice, this direct inversion should be prohibited. Methods such as Gaussian elimination (LU factorization), Gauss-Jordan elimination, Cholesky factorization, or Crout’s method for LU decomposition are available.