|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.
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
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
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];
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.
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.
function F = rope_shoot ( alpha ) % F = rope_shoot ( alpha ) % comments % your name and the dateand this code should do the following:
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 U=inv(A)*b but tells Matlab to solve the equation A*U=b without actually forming the inverse of A.
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.
function [x,U] = rope_bvp(npts) % [x,U] = rope_bvp(npts) % comments % your name and the date
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.
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 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
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(:).
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;
dx=1/(NPTS+1); x=dx*(1:NPTS)'; %column vectorand we will assume an initial temperature distribution
dx=1/(NPTS+1); x=dx*(1:NPTS)'; plot(x, U(k,:) )Please include this plot (for some choice of k>1) with your summary.
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.
plot(x,U(1,:)); hold on for k=2:20 plot(x,U(k,:)); end hold offPlease include this plot with your summary.
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
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.