\magnification=\magstep1
\centerline{Mathematics 174 -- Applied Mathematics 2}
\centerline{Lab/Problem Set 1 -- A ``Shooting Party''}
\centerline{January 27, 2000}
\bigskip
\noindent
{\it Background}
\bigskip
We have now seen that the Sturm-Liouville differential operator
$$L(\phi) = {d\over dx}\left(p(x) {d \phi\over dx}\right) + q(x)\phi,$$
restricted to the space of $C^2$ functions $\phi$ satisfying the 
boundary conditions
$$\cases{\beta_1 \phi(a) + \beta_2 \phi'(a) = 0&\cr
\beta_3 \phi(b) + \beta_4 \phi'(b) = 0&\cr},$$
is {\it self-adjoint\/} for the inner product 
$$\langle f,g\rangle = \int_a^b f(x)g(x)\ dx.$$
Recall that {\it eigenfunctions} are non-zero solutions of
the equation $L(\phi) = -\lambda \sigma(x) \phi$
and the boundary conditions.  

The self-adjointness implies that eigenfunctions
for distinct eigenvalues are orthogonal with respect to the integral
inner product with weight function $\sigma(x)$.  That is, if
$$\cases{L(\phi) = -\lambda \sigma(x)\phi(x) & and \cr
         L(\psi) = -\mu \sigma(x)\psi(x), & where\cr}$$
$\lambda \ne \mu$, then we have
$$\langle\phi,\psi\rangle_\sigma = \int_a^b \phi(x)\psi(x)\sigma(x)\ dx = 0.$$
Also, assuming the Sturm-Liouville problem we consider is
regular, there is a 1-dimensional space of eigenfunctions for 
each eigenvalue.

The other portions of the the Sturm-Liouville theorem (the form of the
set of eigenvalues, the completeness of the eigenfunctions, and 
the pointwise convergence theorem for eigenfunction expansions) are 
more difficult to prove.  We will look at a portion of the proof
that the eigenvalues form a sequence $\lambda_n \to \infty$ in the 
problems below, but we will not go into the details on the remaining
portions.  

Instead, today, we want to illustrate the theorem by looking at an
example where the solutions of the Sturm-Liouville boundary value
problem cannot be expressed in terms of the elementary functions
from calculus.  We will use numerical methods to graph approximate 
solutions and find the eigenvalues.

\vfill\eject
\noindent
{\it Lab Question}
\bigskip
We will consider the Sturm-Liouville problem
$$\cases{{d^2\phi \over dx^2} + \lambda (x^2 + 1)\phi = 0 &\cr
         \phi(0) = 0 &\cr
         \phi(1) = 0 &\cr}$$
where $p(x) = 1$, $q(x) = 0$, and $\sigma(x) = x^2 + 1$.  Since $p(x)$
and $\sigma(x)$ are continuous and strictly positive on $[0,1]$, this is a
regular problem.  Our goals are to approximate the {\it first four eigenvalues}
$\lambda_1,\lambda_2,\lambda_3,$ and $\lambda_4$,
and to generate graphs of corresponding eigenfunctions
$\phi_1,\phi_2,\phi_3,$ and $\phi_4$.  To do this we
will use a rudimentary version of a numerical technique called the 
``shooting method''.  The reason for the name should become clear as
we proceed(!)  

The basic idea of ``shooting'' is that we can solve the differential equation
with {\it initial conditions} $\phi(0) = 0$ and (for instance) $\phi'(0) = 1$
for any given $\lambda$.  Then, by adjusting $\lambda$, we can 
try to get the boundary condition $\phi(1) = 0$ to hold as well.  
(Think of firing a cannon to hit a target at $x = 1$, where you can adjust the
range using the value of $\lambda$!!)

\item{1.} For instance use the following Maple commands to load the {\tt DEtools}
package and define a function called shoot:
\medskip
\centerline{\tt with(DEtools);}
\centerline{\tt shoot := lambda -> DEplot(diff(phi(x),x\$2) + lambda*(x\^{}2+1)*phi(x),}
\centerline{\tt phi(x),x=0..1,[[phi(0)=0,D(phi)(0)=1]],stepsize=0.001);}
\medskip
\item{} Try entering the command {\tt shoot(5)}.  The output is a
plot of the solution of the
differential equation with $\lambda = 5$.  Then enter {\tt shoot(10)}.
What do these graphs tell you about the location of the first eigenvalue:
$\lambda_1$?  Try zeroing in to find a good-looking approximation to 
$\lambda_1$.  

\item{2.} Continue in the same way
to generate approximate values of $\lambda_2,\lambda_3,$
and $\lambda_4$ together with graphs of eigenfunctions for each.  Note
the zeroes of the eigenfunctions in $[0,1]$.

\item{3.}  Of course, the ``eye-ball'' method you used above
is not terribly accurate.  If we wanted better approximations
to the eigenvalues, we would want to be able to examine the actual
approximate value of the computed solution at $x = 1$ and see whether
it is within some specified distance of $0$ (e.g.
$|\phi(1)| < 10^{-d}$ for some $d$).  To do this in
Maple, we can proceed as follows.  First, we can define the 
solution of the initial value problem as a function (depending on 
$\lambda$, and also on $x$):
\medskip
\centerline{\tt sol := lambda -> dsolve($\{$diff(phi(x),x\$2) + lambda*(x\^{}2+1)*phi(x),}
\centerline{\tt phi(0)=0,D(phi)(0)=1$\}$,phi(x),numeric);}
\medskip
\item{}  Then, for instance, {\tt phi := sol(5);} defines the function 
$\phi$ which solves the differential equation and the initial conditions
$\phi(0) = 0$, $\phi'(0) = 1$, with $\lambda = 5$.  If you enter {\tt phi(1);}
the output
$$[x = 1, \phi(x) = .2206662798196227, \phi'(x) = -1.063068688690081]$$
gives the values of $\phi$ and its derivative at $x = 1$.
From this we can see that $\lambda = 5$ is not an eigenvalue!
Take your approximate eigenvalues from parts 1 and 2 and ``polish''
them a bit to get better approximations where the corresponding approximate
eigenfunctions have $|\phi(1)| < 10^{-4}$.
\bigskip
\noindent
{\it Problem Set Question}
\bigskip
In this question, you will work through a big portion of a
proof that a regular Sturm-Liouville 
boundary value problem has a sequence of real eigenvalues $\lambda_n$
with $\lambda_n \to \infty$ as $n \to \infty$.
\bigskip
\item{1)}  (Preparation -- first order linear ODE).  A first order
linear ODE is a differential equation of the form
$${dy\over dx} = a(x) y + b(x)\leqno(1)$$
\itemitem{a.} The corresponding homogeneous equation is 
$${dy\over dx} = a(x) y\leqno(2)$$
show that the general solution of the
homogeneous equation is 
$$y_h(x) = y(0) e^{\int_0^x a(t)\ dt}.$$
\itemitem{b.}  Show that if $y_1(x)$ and $y_2(x)$ are any two solutions
of (1), then $y_1(x) - y_2(x)$ is a solution of the homogeneous
equation (2).
\itemitem{c.}  Show that if we multiply both sides of 
$${dy\over dx} - a(x) y = b(x)$$
by $e^{-\int a(x)\ dx}$, then the resulting equation can be integrated
and solved for $y$ to yield a particular solution of (1):
$$y_p(x) = e^{\int a(x)\ dx} \int e^{-\int a(x)\ dx} b(x)\ dx$$
\itemitem{d.}  From parts a, b, c, deduce that the general solution of
(1) can be written:
$$y(x) = e^{\int_0^x a(t)\ dt}\left(y(0) + \int_0^x e^{-\int_0^s a(t)\ dt} b(s)\ ds\right)$$
\bigskip
\item{2)}  Consider the general (regular) Sturm-Liouville boundary
value problem:
$$\cases{L(\phi) = {d\over dx}\left(p(x) {d \phi\over dx}\right) + q(x)\phi
=-\lambda \sigma(x)\phi & \cr
\beta_1 \phi(a) + \beta_2 \phi'(a) = 0&\cr
\beta_3 \phi(b) + \beta_4 \phi'(b) = 0&\cr}$$
Introduce notation $X(x) = \phi(x)$ and $Y(x) = p(x)\phi'(x)$.  Show that the 2nd order
differential equation can be rewritten as a system of two 1st order
equations as follows:
$$\cases{ X' = {1\over p(x)} Y &\cr
          Y' = -(q(x) + \lambda \sigma(x))X&\cr}$$
\bigskip
\item{3)} Make the following change to ``polar-like'' coordinates:
$X = R \sin\theta, Y = R \cos\theta$. 
(Note: $R$ and $\theta$ are corresponding functions of $x$.)
Show that the system from part 2 becomes
$$\cases{
R' \sin\theta + R \cos\theta\cdot \theta' = {1\over p} R \cos\theta &\cr
R' \cos\theta - R \sin\theta\cdot \theta' = -(q + \lambda\sigma)R \sin\theta &\cr}$$
From this deduce the following two ODE's:
$$\cases{
R' = {R\over 2}\left({1\over p} -(q + \lambda\sigma)\right)\sin 2\theta & (3)\cr
\theta' = {1\over p}\cos^2\theta + (q + \lambda\sigma)\sin^2\theta & (4)\cr}$$
\item{4)}  Equation (4) does not depend on $R$.  Hence general theory for
differential equations implies that for each $\lambda \in {\bf R}$ (and any
fixed $c \in {\bf R}$), there
is a solution $\theta_\lambda(x)$ of (4) with $\theta_\lambda(a) = c$.
Using this, show that the solution of (3) can be written as
$$R(x) = R(a) e^{\int_a^x F_\lambda(t)\sin 2\theta_\lambda(t)\ dt}$$
where $F_\lambda(x) = {1\over p(x)} - (q(x) + \lambda\sigma(x))$.
Deduce that if $R(a) \ne 0$,
then $X(x) = 0$ only when $\theta_\lambda(x) = k\pi$ for some integer $k$.
\item{5)}  So, from part 4, it is enough to show that for each $n \ge 1$, there is
a $\lambda_n$ such that $\theta_{\lambda_n}(b) = n\pi$.  Those give the 
eigenvalues.  (Note that this is essentially what were were doing in 
the lab portion with ``shooting''.)  The key step in showing this is the 
following Lemma:

\bigskip
\proclaim Lemma.  If $\lambda < \mu$, and
$\theta_\lambda(a) \le \theta_\mu(a)$, then $\theta_\lambda(x) < \theta_\mu(x)$
for all $x \in (a,b]$.  

\bigskip
\item{} To show this, let $\omega = \theta_\mu(x) - \theta_\lambda(x)$.
Using equation (4) above (which is satisfied by both $\theta_\mu$
and $\theta_\lambda$), find a 1st order linear ODE
$$\omega'(x) = a(x) \omega(x) + b(x)$$
satisfied by $\omega$, then use part 1 above to determine the form
of the solution with $\omega(a) = \theta_\mu(a) - \theta_\lambda(a) > 0$.
Deduce the Lemma from the form of the solution of your ODE.

\bigskip
\item{} To complete the proof, there are still a few technical details
to handle, but we will not consider them.  The strategy is this:
\itemitem{i.} First, by the Lemma, we know $\theta_\lambda(b)$ is 
increasing as a function of $\lambda$.
\itemitem{ii.}  We would show that $\theta_\lambda(b) < \pi/2$
for some $\lambda$, and that $\theta_\lambda(b) \to \infty$
as $b \to \infty$.  It follows then that for each $n$, there
is some $\lambda_n$ such that $\theta_{\lambda_n}(b) = n\pi$.

\bigskip
\noindent
{\it Assignment}
\bigskip
Individual writeups, due in class on Thursday, February 3.

          	
\bye





