MATH2071: LAB 10: BVPs and PDEs

Introduction Exercise 1
Boundary Value Problems Exercise 2
Shooting Methods Exercise 3
Discretizing a BVP Exercise 4
A Heat Equation Exercise 5
The Method of Lines  


Our study of the initial value problem for ordinary differential equations has taught us how to estimate the behavior, over time, of a scalar quantity, supposing we have an initial value of the quantity, and a description of how it tends to change over time. This very limited ability can nonetheless sometimes be extended to a few other kinds of problems.

One class of problems is known as ``boundary value problems'' (BVPs). A simple example of such a problem would describe the shape of a chain hanging between two posts. We know the position of the endpoints, and we have a second order differential equation describing the shape. If the two conditions were both given at the left endpoint, we'd know what to do right away. But how do we handle this ``slight'' variation?

The other sort of problem that we can take a stab at is the solution of certain ``partial differential equations'' (PDEs), particularly equations that describe the flow of heat. Really interesting problems require more sophisticated treatment, but at least we can get some flavor of this application.

This lab will take two sessions. If you print this lab, you may prefer to use the pdf version.

Boundary Value Problems

A one-dimensional boundary value problem (BVP), is similar to an initial value problem, except that the data we are given isn't conveniently located at a starting point, but rather some is specified at the left end point and some at the right. (We're also usually thinking of the independent variable as representing space, rather than time, in this setting).

Our ideas about solving an initial value problem, (IVP), have to be modified in some way. One possibility is to stubbornly use the IVP idea, and make up any necessary initial conditions. The other idea might be to somehow think of the solution points of the discretized equation as satisfying a set of coupled linear equations that we can solve.

We will be using the following problem as the illustrative example for several of the following exercises.
The clothesline BVP: A rope is stretched between two points. If the rope were weightless, or rigid, it would lie along a straight line; however, the rope has a weight and is elastic, so it sags down slightly from its ideal linear shape. Determine the curve described by the rope. We will use the variable $ x$ to denote horizontal distance and $ u(x)$ to denote height of the rope at the point $ x$.

Forces on a clothesline


The equation for the curve described by the rope can be derived (this is not a proof: it is a description of why you should believe the equation) in the following manner. Suppose that the tension in the rope is $ T$, and consider a tiny piece of the rope of length $ dx$, and with mass per unit length $ {\rho}$. The total mass of the differential piece of rope is $ \rho dx$, so that the force due to gravity is directed downward and is given by This piece of rope observes forces on each of its ends. The magnitudes of these forces are equal to the tension, $ T$, and the directions are given by the slope of the curve at the ends of the differential piece. Hooke's law says that the tension is proportional to the amount of strain in the string, $ T=-Kdu/dx$, where $ K$ is a constant of proportionality. Hence, the equation can be written as

$\displaystyle -K\left.\frac{du}{dx}\right]_{\mbox{right}}+K\left.\frac{du}{dx}\right]_{\mbox{left}}
=- g\rho dx. $

Dividing both sides by $ dx$ and letting $ dx\rightarrow0$ yields the equation

$\displaystyle -\frac{d}{dx}\left(K\frac{du}{dx}\right)=-\rho g$

or, assuming $ K$ is constant in space,

$\displaystyle u''=\rho g/K.$

For the sake of definiteness, take $ g\rho/K=0.2$, the left end of height $ 1$ at $ x=0$, and the right end height of $ 1.5$ at $ x=5$. Thus, the system to be solved is

$\displaystyle u''$ $\displaystyle =$ $\displaystyle 0.2$  
$\displaystyle u(0)$ $\displaystyle =$ $\displaystyle 1$ (1)
$\displaystyle u(5)$ $\displaystyle =$ $\displaystyle 1.5$  

The ODE is linear. Linearity implies good things such as the existence and uniqueness of solutions.

As you know, the ODE (1) can be written as a system of first order differential equations. Consider, as an example, a fourth-order differential equation

$\displaystyle 5u''''+4u'''+3u''+2u'+u=\sin(x)$

with initial conditions

$\displaystyle u_(0)=0,\quad u'(0)=10,\quad u''(0)=20,\quad u'''(0)=30.$

The first step is to define new variables
$\displaystyle y_1$ $\displaystyle =$ $\displaystyle u$  
$\displaystyle y_2$ $\displaystyle =$ $\displaystyle u'$  
$\displaystyle y_3$ $\displaystyle =$ $\displaystyle u''$  
$\displaystyle y_4$ $\displaystyle =$ $\displaystyle u'''$  

(Note that the original equation is fourth order and there are four new variables defined. With these new variables, the original scalar equation can be written as a system of equations
$\displaystyle y'_1$ $\displaystyle =$ $\displaystyle y_2$  
$\displaystyle y'_2$ $\displaystyle =$ $\displaystyle y_3$  
$\displaystyle y'_3$ $\displaystyle =$ $\displaystyle y_4$  
$\displaystyle y'_4$ $\displaystyle =$ $\displaystyle (\sin x-y_1-2y_2-3y_3-4y_4)/5$  

Applying this approach to (1) by identifying $ y_1=u$ and $ y_2=u'$ yields the system

$\displaystyle y'_1$ $\displaystyle =$ $\displaystyle y_2$  
$\displaystyle y'_2$ $\displaystyle =$ $\displaystyle 0.2$  

The right side of this system can be evaluated using the following Matlab code, which should be placed in a Matlab function m-file named rope_ode.m, or downloaded directly as rope_ode.m.
function yprime=rope_ode(x,y)
% yprime=rope_ode(x,y) computes the
% rhs of the first-order system corresponding to u'' = 0.2

% your name and the date
yprime=[ y(2)

Shooting Methods

Since we don't have any idea how to solve a BVP, but we know how to solve an IVP, let's just do what we do best. Taking as our example the rope BVP, how much violence do we have to do in order to make it look like an IVP? Well, we expect to have two conditions at the left point and none at the right point. So let's temporarily consider the related problem, where we have made up an extra boundary condition at our initial value of $ x=0$:

$\displaystyle u''&=0.2$    
$\displaystyle u(0)$ $\displaystyle =1$ (2)
$\displaystyle u'(0)$ $\displaystyle =\alpha$    

The following pair of exercises attacks this system using a method called ``shooting.'' The strategy behind this method is the following.

For this pair of exercises, you will be using built-in Matlab ODE routines to solve the systems of ODEs that arise.

Exercise 1:
  1. Copy the above code into a file named rope_ode.m. Note that since we plan to use ode45 there is no need to add the Jacobian matrix to rope_ode.m! Matlab computes this Jacobian internally. This is a great convenience, but in many cases you would have to provide a function for the Jacobian or ode45 (or ode15s, etc.) might fail. In that case, setting an option allows the Jacobian computation.
  2. Choose a provisional value $ \alpha=0$ and use the Matlab function ode45 to solve the system (2) on the interval [0,5]. What is the value of $ u(5)-1.5$? Is it positive or negative? (Be careful: which of the components of y coesponds with $ u$?)
  3. By trial and error, find some value of $ \alpha$ for which the value of $ u(5)-1.5$ is of the opposite sign as for $ \alpha=0$.

Any nonlinear equation solver can be used to systematically solve for the value of $ {\alpha}$ so that $ F(\alpha) = 0$, and each has its idiosyncrasies. You have just found two values of $ \alpha$ for which the right-endpoint error has opposite signs, so something like bisection would work nicely. Matlab already provides a function, fzero that will find the solution of a function that changes sign, and you can use that one.

Exercise 2:
  1. Write an m-file called rope_shoot.m that accepts a value of alpha and evaluates F(alpha), for the rope BVP. The file should have the signature
    function F = rope_shoot ( alpha )
    % F = rope_shoot ( alpha )
    % comments
    % your name and the date
    and this code should do the following:
    • Use the input value of alpha as the initial condition for $ y'(0)$;
    • Use ode45 to compute the solution [x,y] of the IVP (2) defined by the initial conditions, and the right hand side function rope_ode, for $ x\in[0,5]$;
    • Return in the function value F the value of $ u(5)-1.5$. If you call the Matlab results from ode45 x and y, then F=y(end,1)-1.5. (Recall that end is a special Matlab variable that means the largest value of the subscript.)
    For a given value of $ \alpha$, the function you just wrote will return $ y(5)-1.5$. When $ \alpha$ is just right, it will return 0.
  2. Use the Matlab function fzero to find the value of $ \alpha$ that makes $ F(\alpha)=y(5)-1.5=0.$ fzero requires two parameters, the name of a function first and second the vector [alpha1,alpha2] of the two values of $ \alpha$ that you found in Exercise 1 for which $ F$ has opposite signes. What is the value of alpha you found?
  3. Plot the solution you found. Does the curve have a height of 1 at $ x=0$ and a height of 1.5 at $ x=5$? You do not need to send me this plot.

Discretizing a BVP

The ``derivation'' presented above for the shape of the rope is suggestive of another way to solve for the shape. Assume that we have divided the interval up into $ N$ equal intervals of width $ \Delta x$ determined by $ N+1$ points. Denote the spatial points $ x_{n}$, $ n=0,1,\ldots,N+1$. Approximate the value of $ y(x_{n})$ by $ y_{n}$. Now, consider the $ n^{\mbox{th}}$ interval as if it were the differential piece of rope mentioned in the derivation. Using the standard approximation for a derivative, the slope of the rope at the left of the $ n^{\mbox{th}}$ interval could be approximated as $ (y_{n}-y_{n-1})/\Delta x$ and the slope on the right of the $ n^{\mbox{th}}$ interval could be approximated as $ (y_{n+1}-y_{n})/\Delta x$. Put all these together and you get

$\displaystyle \frac{y_{n+1} - 2 y_{n}+ y_{n-1}}{\Delta x^{2}} = \frac{g\rho}{K} = 0.2

We can associate this equation with the solution value at $ y_{n}$, except for $ n=0$ and $ n=N+1$ (do you see why). Conveniently, those are the points at which we have boundary conditions specified.

In particular, let us look at approximating our rope BVP at 6 points. We set up the ODE at points 1, 2, 3, and 4, and associate the boundary conditions with the $ n=0$ and $ n=5$ solution values. I also multiplied through by the denominator to make things look nicer:

$\displaystyle u_0$ $\displaystyle =$ $\displaystyle 1$  
$\displaystyle u_{0} - 2 u_{1} + u_{2}$ $\displaystyle =$ $\displaystyle 0.2 \Delta x^{2}$  
$\displaystyle u_{1} - 2 u_{2} + u_{3}$ $\displaystyle =$ $\displaystyle 0.2 \Delta x^{2}$ (3)
$\displaystyle u_{2} - 2 u_{3} + u_{4}$ $\displaystyle =$ $\displaystyle 0.2 \Delta x^{2}$  
$\displaystyle u_{3} - 2 u_{4} + u_{5}$ $\displaystyle =$ $\displaystyle 0.2 \Delta x^{2}$  
$\displaystyle u_{5}$ $\displaystyle =$ $\displaystyle 1.5$  

Actually, in Equation (3), the quantities $ u_0$ and $ u_5$ are not really variables, being fixed by the boundary conditions. Hence the only variables are $ u_1$, $ u_2$, $ u_3$ and $ u_4$. The system can be rewritten as

-2u_1&+u_2&+0&+0&=&.1\Delta x^2-u_0 ...
...&.1\Delta x^2\\
0&+0&+u_3&-2u_4&=&.1\Delta x^2-u_5

and this system has been formatted to suggest the matrix equation

$\displaystyle \left[\begin{array}{rrrr} -2&1&0&0\\ 1&-2&1&0\\ 0&1&-2&1\\ 0&0&1&...
...elta x^2-u_0\\ .1\Delta x^2\\ .1\Delta x^2\\ .1\Delta x^2-u_5\end{array}\right]$ (4)

By discretizing the differential equations we have created a set of linear alebraic equations that have the symbolic form $ AU=b$. To set up and solve the equations (4) in Matlab, we could type:

NPTS = 4;
A = [ -2  1  0  0;
       1 -2  1  0;
       0  1 -2  1;
       0  0  1 -2];
% NPTS interior mesh points, NPTS+1 intervals
dx = 5.0 / ( NPTS + 1 );
x = dx * (1:NPTS);
b = [ 0.2*dx^2-uLeft
U = A \ b
Make sure you understand the first and last components in b. You should recall that the backslash notation is shorthand for saying U=inv(A)*b but tells Matlab to solve the equation A*U=b without actually forming the inverse of A.

Exercise 3: In this exercise, you will be using the above code to solve the rope BVP.
  1. Copy the above code to a script m-file named ex3.m, and execute it to find a solution of Equation (4).
  2. Verify by hand that the values of U and b that you found satisfy at least one of the linear equations. You do not have to send me your handwritten verification.
  3. Then plot U versus x. Include the printed values of U as part of the lab summary, but you don't have to send the plot. When you examine the plot, don't forget that the end points, x=0 and x=5 are not included.

The purpose of the previous exercise is to verify that you have copied the code correctly. You should always use a small, simple problem to verify code by comparison with hand calculations. In the next exercise you will be modifying the above code to handle the case of large NPTS and solving a slightly more realistic problem. Of course, you don't want to bother typing in the matrix A if NPTS is 100 or more, so you will be writing Matlab code to do it.

Exercise 4:
  1. Modify the script m-file ex3.m into a function m-file called rope_bvp.m with the signature
    function [x,U] = rope_bvp(npts)
    % [x,U] = rope_bvp(npts)
    % comments
    % your name and the date
  2. Add comments after the signature line and modify the matrix (A) and right side vector (b) generation statements to be valid for arbitrary values of npts. (You should also eliminate the line NPTS = 4.) Hint: You can use the zeros(npts,npts) statement to generate an npts-by-npts matrix of all zeros for A and then fill in the non-zero values. You can use the ones(npts,1) to construct a column vector of length npts containing all ones.

    Got a minute for a puzzle? The function ones(m,n) returns an m by n matrix of 1's. The function diag(v,d) returns a matrix whose d-th diagonal contains the values in the vector v. Using these ideas, it is possible, but not required, to compute the matrix A in a single (logical) line.

  3. Check your work by running rope_bvp for npts=4 and confirming that you get the same values of U as in the previous exercise. Name these variables x1 and U1.
  4. Now we are ready to solve the big problem! Use npts=100 so that there are 100 unknowns U(1:100). Call the new solution [x2,U2], and plot U2 versus x2. Plot U1 versus x1 as circles (plot(x1,U1,'o') on the same frame (hold on) and send the single frame with both plots to me as part of the summary. Part of the point of this exercise is for you to see the effect of increasing accuracy in this kind of problem.
  5. Use a finite difference expression for the derivative to estimate the derivative of U at the left endpoint (don't forget that the value of the solution at the endpoint itself is the boundary value and that U(1) is for x=dx). How does it compare with the value you computed in Exercise 2?

Remark: You may wonder why the four-point mesh solution and the hundred-point mesh solution agree at the four common points. This behavior is highly unusual. It happens because the exact solution of the differential equation is a quadratic and the difference scheme exactly reproduces quadratic functions.

A Heat Equation

A partial differential equation (PDE) involves derivatives of a function $ u$ which depends on more than one independent variable. One of the classic PDEs is the one-dimensional heat equation, which describes the behavior of many physical quantities, including the temperature $ u(x,t)$ of a rod. You could imagine a rod of some sort of metal, part of which might be heated in some manner (held over a flame for instance), whose ends are kept at a constant temperature (encased in ice for instance), and which starts out with some distribution of temperature. The partial differential equation describing the temperature $ u(x,t)$ as a function of time is

$\displaystyle \frac{\partial u}{\partial t} = \frac{\partial^{2}u}{\partial x^{2}} + f(x,t)

where the function $ f(x,t)$ describes the heat generation, which might change with time and space. Boundary conditions describing the effect of the ice (0 degrees) at the end points can be written as
$\displaystyle u(0,t)$ $\displaystyle =$ 0  
$\displaystyle u(1,t)$ $\displaystyle =$ 0  

and we also specify an initial condition, that is, the temperature of every point on the rod at the starting time,

$\displaystyle u(x,0) = g(x)

Of course, other boundary and initial conditions are possible.

The Method of Lines

Look at the system and pretend, just for a moment, that time is frozen. The system suddenly starts looking just like one of the boundary value problems we just solved! This observation is the basis for the ``method of lines,'' wherein the spatial discretization is performed separately from the temporal. Consider the function $ u(t,x)$, but think of it as a function of $ x$ first. The resulting BVP (ignoring the heat generation rate $ f(x,t)$ for simplicity) comes from the heat equation, written with primes denoting spatial differentiation so it looks more like what we have been doing.

$\displaystyle u''=\frac{\partial u}{\partial t}$

along with boundary conditions at, say, $ x=0$ and $ x=1$:

$\displaystyle u(0,t)=u(1,t)=0.$

We solved a BVP a lot like this one above. We broke the interval $ [0,1]$ into $ {\tt npts}+1$ subintervals and labelled the $ ({\tt npts}+2)$ resulting points $ x_{n}, n=0,1,2,...,({\tt npts}+1)$. Next, we defined $ u_{n}$ as being the approximate solution at $ x_{n}$. Keep in the back of you mind, though, that $ u_{n}$ is really still a function of $ t$, $ u_{n}(t)$. We will denote the vector of values $ u_{n}$ by $ U$, because we have already used the unsubscripted $ u$ to denote the continuous solution. Keep in mind that $ U(t)$ is a function of time. We constructed the tridiagonal matrix $ A$, and the resulting system of equations becomes

$\displaystyle \frac{dU}{dt}=(\frac1{\Delta x^2})AU.$

But we remain on familiar ground: this is just a system of IVPs! We know a bunch of different methods to solve it. It turns out that the system is moderately stiff, and will become more stiff when npts is taken larger and larger. We will use backwards Euler to solve this system. Recall that backwards Euler requires an m-file to evaluate the partial derivative as well as an m-file to evaluate the function.

In Matlab notation, the variable $ U$ will be a matrix whose entries $ U_{kn}$ approximate the values $ u_n(t_k)$. In Matlab notation, for the $ k^{\mbox{th}}$ time interval, U(k,:) represents the (row vector of) values at the locations x(:). The initial condition for U should be a column vector whose values are specified at the locations x(:).

Exercise 5:
  1. Write a function m-file called heat_ode.m to construct the spatial discretization of the heat equation and its Jacobian. This discretization is the same as rope_bvp.m, but the zero boundary conditions make it even simpler. Here is an outline:
    function [Udot,J]=heat_ode(t,U)
    % [Udot,J]=heat_ode(t,U)
    % compute the right side of the time-dependent ODE arising from
    % a method of lines reduction of the heat equation
    % and its Jacobian derivative, J.
    % Boundary conditions are fixed =0 at the
    % endpoints x=0 and x=1.
    % A fixed number of spatial points (=NPTS) is used.
    % The variable t is not used, but is kept as a place holder.
    % U is the vector of the approximate solution at all spatial points
    % output Udot is the time derivative of U
    % output J is the Jacobian matrix of
    %   partial derivatives of Udot with respect to U
    % spatial intervals
    % Construct the NPTS-by-NPTS A-matrix, just like rope_bvp.m
    A = ???
    % dU/dt is returned in the variable Udot
    Udot = (???)*A*U;
    % The Jacobian (d Udot/dU) is
    J = (???)*A;
  2. The spatial values represent a uniform mesh
    x=dx*(1:NPTS)';   %column vector
    and we will assume an initial temperature distribution
  3. Retrieve the copy of back_euler.m you wrote for last lab or download my copy.
  4. Use back_euler to solve the heat equation starting from UInit You should use 60 steps from time t=0 to time t=1.5. What is the value of U(10,20)?
  5. Recall that the rows of U represent temperatures at different places all at the same time and columns of U represent temperatures at different times all at the same place. To see a snapshot of the solution at timestep k, use the command
    plot(x, U(k,:) )
    Please include this plot (for some choice of k>1) with your summary.
  6. To see the time evolution of the solution value at a fixed point n, use the command
    plot( t,U(:,n) )
    The time evolution is almost exponential decay, once it gets going. You can see this with the command (for, say n=10)
    plot( t,log10(U(:,n)) )
    Please include these plots with your summary.
  7. For this problem, the solution is pretty simple, and you can see its evolution nicely using the following steps
    hold on
    for k=2:20
    hold off
    Please include this plot with your summary.
  8. You can also see a ``flicker picture'' of the evolution with the following steps
    for k=2:60
  9. One way to get a nice plot of the entire solution is to type
    mesh (x,t,U )
    and if you want to highlight the approximately exponential decay of this solution in time, you can use
    mesh (x,t,log10(U) )
    Send me copies of your plots as part of your lab summary.

Back to MATH2071 page.

Mike Sussman 2009-02-04