Classroom Tips and Techniques:
Numeric Solution of a TwoPoint BVP
Robert
J. Lopez
Emeritus Professor of Mathematics and Maple Fellow
Maplesoft
Introduction
A recent query from a colleague in engineering posed a nonlinear twopoint boundary value problem that could only be solved numerically. In this month's article, we'll explore the solution of a similar problem, highlighting the way Maple can be used to obtain the required solution.
Initializations
>


>


Statement
of the BVP
Consider
the nonlinear ODE
>


and the related nonlinear boundary conditions
>


that pose a twopoint boundary value problem on the interval
Using dsolve
to Obtain a Solution
A naive
use of Maple's dsolve
command did not provide a solution to the problem originally posed, running
at great length and then crashing. However, inclusion of the approxsoln
parameter very quickly led to the desired solution. This parameter provides
an initial guess for an iterative RayleighRitz solution process that converges
very quickly to a solution. In the original engineering problem there
were two possible solutions, only one of which was physically meaningful.
If we
apply that same strategy to the problem posed in the previous section, assuming
that perhaps a solution might be close to ,
we obtain the graph shown in Figure 1.
>



Figure
1
Solution of nonlinear twopoint BVP starting with initial
guess

In principle, it should be possible to solve this BVP via a shooting method,
using the two boundary conditions as the governing equations that determine
and .
In the sequel, we apply this idea to the solution process, and also
come to grips with another disturbing facet of the problem we selected, namely,
finding all possible solutions.
A
Solution from First Principles
We begin
by defining a procedure
that, given values for ,
returns the list
This procedure implements a shooting method, starting from the initial
values of
and ,
numerically integrating a solution for the resulting initial value problem,
and returning the endpoint values.
>

h := proc(a,b)
local F;
if type(a, numeric) or type(b, numeric) then
F := dsolve({q, y(0)=a,D(y)(0)=b}, y(x), numeric);
[op([2,2],F(1)),op([3,2],F(1))];
else 'h'(a,b);
end if;
end proc: 
In terms of and
,
the boundary conditions can be expressed via the expressions
>


It is now possible to graph the functions
and ,
resulting in Figure 2.
>



Figure
2
Graph of the equations for
and ,
as determined by the boundary conditions of the BVP

Figure 2 suggests that there are actually four solutions of the given BVP.
Unfortunately,
the fsolve command in Maple 11 does not permit systems of procedure calls,
as would be required by our definitions of
and
through the procedure
Fortunately, this shortcoming is corrected in the next version of Maple.
Meanwhile, we need to solve a pair of nonlinear algebraic equations
whose terms are calculated via the solution of an initial value problem.
We propose
three methods for this. First, we implement Broyden's method as a generalization
of the Secant algorithm for solving a single equation. Second, we apply
a minimization technique to ,
a sum of squares. Finally, because of the special nature of the first
boundary condition, we can eliminate
from the second equation, and apply the Secant method to the single resulting
equation.
In anticipation
of all three iterative techniques, we define the following initial points,
each taken from Figure 2.
>


Broyden's
Method
Newton's
method for solving the system of equations is
captured by the fixedpoint scheme
where
is the Jacobian matrix for F.
Since this requires computation of derivatives, it would be tedious
to implement a method for repeatedly computing .
As an alternative, we consider Broyden's method wherein
is only computed once. Thereafter, an approximation to
is computed, an approximation that converges to where
is a solution to the system
After
obtaining
either analytically or approximately by numeric differentiation, Broyden's
method thereafter replaces with
as an
approximation to
Our implementation
of Broyden's method takes the vector
as impute, along with an initial point, a stepsize
used for numerically computing ,
and the integer
used to define the computational tolerance .
Our code is simplistic, in that it does not typechecking, and does
not make provision for generalization to other systems. It is provided
as a sketch to show how such a computation might be implemented.
>

Broyden
:= proc(F,X0,d,N)
local X, F0, F1, f, g, fa, ga, fb, gb, A0, A1, x0, x1, tol, k, L, Z;
x0 := X0;
X := <a,b>:
f := F[1];
g := F[2];
Z := Equate(X,x0);
F0 := eval(F,Equate(X,x0)):
fa := (eval(f,Equate(X,x0+<d,0>))  F0[1])/d:
ga := (eval(g,Equate(X,x0+<d,0>))  F0[2])/d:
fb := (eval(f,Equate(X,x0+<0,d>))  F0[1])/d:
gb := (eval(g,Equate(X,x0+<0,d>))  F0[2])/d:
A0 := Matrix([[fa,fb],[ga,gb]]):
x1 := x0  A0^(1).F0:
Digits := 30:
tol := 10^(N):
for k from 1 to 10 do
F1 := eval(F, Equate(X,x1));
L := x1  x0;
A1 := A0 + (F1  F0  A0.L).L^%T/(L.L);
Z := A1^(1).F1;
if sqrt(Z.Z)<tol then break;
else x0 := x1; x1 := x1Z; A0 := A1; F0 := F1;
end if;
end do:
Digits := 10:
x1;
end proc: 
Using our Broyden procedure, we compute solutions for
and .
>


With pairs
determined, we can solve the corresponding four initial value problems, and
graph their solutions. The following commands generate the solutions
and graphs,
>


and the graphs themselves appear in Figure 3.
>



>



>



>



Figure
3
Graphs of four possible solutions generated as solutions
of initial value problems from pairs
computed by Broyden's method

The four solutions found in Figure 3 can be obtained directly with Maple's
dsolve command, provided an appropriate approximate solution is provided.
If we list the four approximate solutions
>


then the following loop computes the corresponding dsolve
solutions.
>


Figure 4 displays these solutions.
>



>



>



>



Figure
4
Solutions obtained by dsolve
using the approxsoln parameter

The graphs in Figures 3 and 4 are essentially the same.
Minimizing
a Sum of Squares
A pair
that minimizes the sum of squares
>


is a solution to the equations .
We can compute such minimizing pairs with the NLPSolve
command in Maple's Optimization
package. This command provides the nonlinear simplex (NelderMead) method,
which does not require computation of derivatives. The computed pairs
and the values of
at the minima are given by the following loop.
>


The first number returned in each solution is the value of
at the computed extremum.
To facilitate
a comparison with the values computed by Broyden's method, we express the
nonlinear simplex solutions as the vectors
>


The Euclidean norms of the differences between the two solutions are given
by
>


The agreement is substantial.
Secant
Method
If the
first boundary condition is solve for ,
we obtain
>


with which
can be eliminated from the second boundary condition. Thus, we have
the single equation
>


to solve. Figure 5 provides a graph of the lefthand side of this equation.
>



Figure
5
Graph of

There are four intercepts for the graph in Figure 5, consistent with our earlier
findings on the number of solutions of the nonlinear BVP. The equation
can be solved with the Secant method, which we have coded as the following
procedure. Our procedure includes a calculation of ,
and returns both
and .
>

Secant
:= proc(X0,N)
local y0,k,x0,x1,x2,y1,X,Y,tol;
Digits := 30:
tol := 10^(N):
x0 := X0;
y0 := eval(Ga,a=x0):
x1 := 1.05*x0:
unassign('X','Y');
for k from 1 to 10 do
y1 := eval(Ga,a=x1);
x2 := x1  y1 * (x1  x0)/(y1  y0);
if abs(x2x1)<tol then X:=evalf[N](x2); Y:=evalf[N](eval(B,a=x2));break:
else x0:=x1; y0:=y1; x1:=x2; y1:=eval(B,a=x1);
end if;
end do:
Digits := 10:
<X,Y>;
end proc: 
An application of the Secant method leads to the following four solutions.
>


As we did earlier, we compare these solutions with those computed by Broyden's
method.
>


The solutions
appear to be identical.
Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for noncommercial, nonprofit use only. Contact the author for permission if you wish to use this application in forprofit activities.