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 |

Introduction

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 to denote
horizontal distance and to denote height of the rope
at the point .

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 , and consider a tiny piece of the rope of length , and with mass per unit length . The total mass of the differential piece of rope is , 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, , 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, , where is a constant of proportionality. Hence, the equation can be written as

For the sake of definiteness, take
, the left end of
height at , and the right end height of at .
Thus, the system to be solved is

The ODE is

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

(

Applying this approach to (1) by
identifying and yields the system

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

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) 0.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 :

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

- For each value of , we can solve this problem, and get a numerical solution at a sequence of points up to the right endpoint.
- Since every value of determines a (numerical) solution , we can regard the difference between the value we got, and the value we want, as a function . (Where we write , but we should really write something like , to emphasize that the solution depends on the parameter.)
- The BVP solution we are looking for has the property that .

**Exercise 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. - Choose a provisional value and use the Matlab
function
`ode45`to solve the system (2) on the interval [0,5]. What is the value of ? Is it positive or negative? (Be careful: which of the components of`y`coesponds with ?) - By trial and error, find some value of for which the value of is of the opposite sign as for .

- Copy the above code into a file named

Any nonlinear equation solver can be used to systematically
solve for the value of so that
, and each has its
idiosyncrasies. You have just found two values of
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**:- 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 signaturefunction 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 ; - 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 ; - Return in the function value
`F`the value of . 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.)

- Use the input value of
- Use the Matlab function
`fzero`to find the value of that makes`fzero`requires two parameters, the name of a function first and second the vector`[alpha1,alpha2]`of the two values of that you found in Exercise 1 for which has opposite signes. What is the value of`alpha`you found? - Plot the solution you found. Does the curve have a height of 1 at and a height of 1.5 at ? You do not need to send me this plot.

- Write an m-file called

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 equal intervals of width determined by points. Denote the spatial points , . Approximate the value of by . Now, consider the 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 interval could be approximated as and the slope on the right of the interval could be approximated as . Put all these together and you get

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 and solution values. I also
multiplied through by the denominator to make things look nicer:

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

By discretizing the differential equations we have created a set of linear alebraic equations that have the symbolic form . 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); uLeft=1; uRight=1.5; b = [ 0.2*dx^2-uLeft 0.2*dx^2 0.2*dx^2 0.2*dx^2-uRight] U = A \ bMake sure you understand the first and last components in b. You should recall that the backslash notation is shorthand for saying

**Exercise 3**: In this exercise, you will be using the above code to solve the rope BVP.- Copy the above code to a script m-file named
`ex3.m`, and execute it to find a solution of Equation (4). - 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. - 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.

- Copy the above code to a script m-file named

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**:- Modify the script m-file
`ex3.m`into a function m-file called`rope_bvp.m`with the signaturefunction [x,U] = rope_bvp(npts) % [x,U] = rope_bvp(npts) % comments % your name and the date

- 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. - 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`. - 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. - 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?

- Modify the script m-file

**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 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
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 as a function of time is

0 | |||

0 |

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

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 , but think of it as a function of first. The resulting BVP (ignoring the heat generation rate for simplicity) comes from the heat equation, written with primes denoting spatial differentiation so it looks more like what we have been doing.

We solved a BVP a lot like this one above. We broke the interval into subintervals and labelled the resulting points . Next, we defined as being the approximate solution at . Keep in the back of you mind, though, that is really still a function of , . We will denote the vector of values by , because we have already used the unsubscripted to denote the continuous solution. Keep in mind that is a function of time. We constructed the tridiagonal matrix , and the resulting system of equations becomes

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 will be a matrix whose
entries approximate the values . In Matlab
notation, for the
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**:- 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 NPTS=100; dx=1/(NPTS+1); % 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;

- The spatial values represent a uniform mesh
dx=1/(NPTS+1); x=dx*(1:NPTS)'; %column vector

and we will assume an initial temperature distributionUInit=x.^2;

- Retrieve the copy of
`back_euler.m`you wrote for last lab or download my copy. - 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)`? - 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 commanddx=1/(NPTS+1); x=dx*(1:NPTS)'; plot(x, U(k,:) )

Please include this plot (for some choice of`k>1`) with your summary. - To see the time evolution of the solution value at a fixed point
`n`, use the commandplot( 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. - For this problem, the solution is pretty simple, and you
can see its evolution nicely using the following steps
plot(x,U(1,:)); hold on for k=2:20 plot(x,U(k,:)); end hold off

Please include this plot with your summary. - You can also see a ``flicker picture'' of the evolution
with the following steps
plot(x,U(1,:)) axis([0,1,0,1]) for k=2:60 pause(0.1); plot(x,U(k,:)); axis([0,1,-.1,1]) end

- 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 usemesh (x,t,log10(U) )

Send me copies of your plots as part of your lab summary.

- Write a function m-file called

Back to MATH2071 page.

Mike Sussman 2009-02-04