Last modified 5 April 2004 by jdavis@math.wisc.edu
In this lab we use MATLAB to implement numerical methods that compute derivatives and solve ordinary differential equations. It is assumed that you remember everything from the first lab; refer back as needed. Please hand in the marked exercises on paper, all together.
Before we get started, we need to know just a bit more MATLAB. The feval function takes in the name of a function, along with some arguments for that function, and evaluates the function with those arguments. For example, the following two commands are equivalent:
sin(0)
feval('sin', 0)
You can also use feval on functions that you yourself write in M-files. For example, suppose you write a function called "add", that takes two arguments and simply adds them. Then the following two commands are equivalent:
add(3, 5)
feval('add', 3, 5)
You might be wondering what the point of all of this is. In later sections, the feval function will be useful in writing MATLAB functions that take other functions as input. For a taste, here's a function that takes in any function and evaluates that function at 0.
function Y = evaluateatzero(F) % takes the name any function F : R -> R as input and returns F(0) Y = feval(F, 0);
Now these commands are equivalent:
sin(0)
evaluateatzero('sin')
feval('evaluateatzero', 'sin')
Suppose we are given a function F : R -> R, written X = F(T), and we want to find its derivative at some fixed T0. The derivative is the slope of the tangent line at T0. We can numerically approximate the derivative by finding the slope of a nearby secant line. That is, we choose some small step size S, and then compute the slope of the line from the point (T0, F(T0)) to the point (T0 + S, F(T0 + S)). This makes for a good warm-up exercise.
Exercise A: Write a function that takes as input (the name of) a function F, a number T0, and a step size S, and returns the approximate derivative of F at T0, as described above. The first line of your function should befunction D = derivative(F, T0, S)
You could make this more accurate (possibly) by averaging the slopes of two secant lines, one through (T0, F(T0)) and (T0 + S, F(T0 + S)), the other through (T0, F(T0)) and (T0 - S, F(T0 - S)). But it seems that most people don't do this.
In the previous section we were given a function and we numerically approximated its derivative. In this section we do the opposite: Given information about a function's derivative, we approximate the function. In particular, we are concerned with differential equations of the form
dX/dT = M(T, X),
as discussed recently in class. Given the function M and an initial condition (T0, X0), we want to find X as a function of T.
Recall Euler's method. First we evaluate M at (T0, X0) to get a slope M0 at that point. Then we travel from (T0, X0), along a line of slope M0, for a certain step size S. We arrive at a new point (T1, X1) = (T0 + S, X0 + S * M0). Starting from this point, we again compute a slope and take a step of size S. We keep repeating this until we have gone as far as we want.
Exercise B: Write a function that takes as input (the name of) a function M, an initial condition (T0, X0), an end value T1, and a step size S. It uses Euler's method to approximate X(T) on the interval from T0 to T1. It returns a 2 x ? matrix, TX, containing the coordinates of all of the step points. That is, TX is a list of points, in that each column in TX is a point. The first point in TX should be (T0, X0). The first line of your function should befunction TX = eulersmethod(M, T0, X0, T1, S)
Now that we have a function that implements Euler's method, we can use it in a variety of ways. For starters, we can visualize how solutions change when we vary the initial condition.
Exercise C: Write a function that takes in the same inputs as the eulersmethod function, except that now the numbers T0 and X0 have been replaced with vectors T and X, giving a list of different initial conditions. For each initial condition in the list, use your eulersmethod function to approximate the solution. Then graph the solutions for all of the initial conditions, together on one big plot. (To accomplish this, use the close command to close the figure window, then hold on, then draw the graphs, then hold off.) Your function should look likefunction ploteulersmethod(M, T, X, T1, S)
Try out this function with M(T, X) = T + X, step size 0.01, t1 = 1, and initial conditions
(-3, 1.5), (-3, 1.7), (-3, 1.9), (-3, 2.1), (-3, 2.3).
Print out your figure (or make a detailed sketch) to hand in, along with your function.
Suppose that we have a grain or body of rock, and that X is the (average) concentration of some mineral/element/isotope in that body. The concentration is decreasing, due to diffusion. For theoretical reasons, we believe that the diffusion follows the differential equation
dX/dT = -X2 + e-T.
This is not realistic; a true diffusion problem is a PDE. Anyway, we have here a nonlinear first-order ODE, not tractable with our earlier methods, but ripe for solving with Euler's method.
Suppose that we know what the concentration is at time T = 1, and we want to know what the concentration was at time T = 0. Well, actually, due to imprecision in our measurements, all we know is that the concentration is somewhere between 0 and 1 when T = 1. The question is: What was the possible range of concentrations at time T = 0?
Exercise D: In an M-file, define the appropriate function M for this ODE. Then, as in Exercise C, try out various initial conditions at time T = 0, until you find the range of initial conditions that produce concentrations between 0 and 1 at time T = 1.Submit the M-file for your function M, along with your answer to the problem.