This assignment was given on 17 March 2009. The work is due 24 March 2009.
Question 1: ARTILLERY, an initial value problem.
The X and Y coordinates of a cannonball are functions of time. We write the position as (X(T), Y(T)). We suppose the cannonball starts out at time T = 0 with (X(0),Y(0)) = (0,0). The cannonball is given an initial velocity V with components (VX,VY). We are measuring in units of seconds and meters, so the gravitational acceleration is 9.8 meters per second squared.
The position of the cannonball is described by the following differential equation, which is called an initial value problem:
X(0) = 0 <---- Initial conditions.
Y(0) = 0
dX/dt = VX <---- Differential equations describe
dY/dt = VY - 9.8 * t how the initial condition changes.
Assume that VX = 10 m/s and VY = 49 m/s. Using calculus, it is possible to write simple exact formulas for the X and Y coordinates of the cannonball.
A) What are these exact formulas?
X(t) = ______________________
Y(t) = ______________________
Using the formula, create a table of the value of Y(t) at one second intervals from t = 0 to t = 10. We will need this later.
If we didn't know calculus, we can't solve the initial value problem exactly. But we can reason as follows. The derivative dY/dt tells us the instantaneous rate of change of Y. If we know the value Y(t), and we want to estimate Y(t+DT), its value DT seconds later, we estimate it using a finite difference approximation
Y(t+DT) = (approximately) Y(t) + Y'(t) * DT
For our data, if we use a time stepsize of DT = 1 second, then our estimate for Y(1) would be:
Y(0+1) = Y(0) + Y'(0) * 1 = 0 + ( 49 - 9.8 * 0 ) * 1 = 49
Now that we have estimate Y(1), we can use the same trick to estimate Y(2):
Y(1+1) = Y(1) + Y'(1) * 1 = 49 + ( 49 - 9.8 * 1 ) * 1 = 88.2
Using a stepsize of DT = 1 second, compute the approximate values of Y for t = 0 up to 10 seconds. Make a table that compares the exact values of Y from the formula, and your finite difference estimate with DT = 1.
Recompute your finite difference estimate by using a stepsize DT = 1/2 second. Compare this estimate of the solution with your exact formula and your previous estimate, using a table or a plot.
B: Turn in a table which lists the values of Y(t) for t = 0, 1, ... to 10, for your exact formula, your DT=1 estimate, and your DT=1/2 estimate.
Question 2: PREDATOR/PREY, a pair of initial value problems.
In this problem, we have two variables and two equations, but now the variables interact, and we must solve the equations together. In such a system, we say the variables are coupled...it just means that changes in one variable affect the other variable.
We suppose that we have a population of harmless bunnies. If left alone, then roughly speaking, each year there would be two more bunnies for each current bunny. If this was all we knew about bunny biology, the equation would indicate exponential growth:
dB/dt = + 2 B(t)
However, there is also a population of foxes F(t). The chances that a bunny will
meet a fox are proportional to B(t)*F(t). We will estimate the effect
of these encounters as a negative factor for the bunny population growth rate:
dB/dt = + 2 B(t) - 0.001 B(t) F(t)
The fox population has two influences on its growth; the fox population tends to starve to death at a constant rate (a negative) but each time a fox encounters a bunny this is GOOD for the foxes (a positive). The starvation is proportional to the number of foxes; the feeding encounters are proportional to the product of bunnies and foxes:
dF/dt = - 10 F(t) + 0.002 B(t) F(t)
Suppose that we know that on January 1, 2010, there are 5000 bunnies and 100 foxes. How will the population evolve over 5 years?
We start by choosing a number of time steps STEP_NUM. Then the size of our time step will be DT = ( 5 / STEP_NUM ), since we're measuring in years. We will try to set up tables T(I), B(I), and F(I), that record, for I from 0 to STEP_NUM, the time, and the number of bunnies and foxes.
Again, the initial conditions help us get started:
T(0) = 0.0
B(0) = 5000
F(0) = 100
(MATLAB won't actually let you use "0" as the index of an array. There are
simple ways to deal with this fact, such as counting from 1 to STEP_NUM+1 instead of
from 0 to STEP_NUM.)
To get the values at the next time step, we use the finite difference idea again:
T(I+1) = T(I) + DT
B(I+1) = B(I) + DT * ( 2 * B(I) - 0.001 * B(I) * F(I) )
F(I+1) = F(I) + DT * ( - 10 * F(I) + 0.002 * B(I) * F(I) )
Write a MATLAB program which takes as input a value for STEP_NUM, and which solves these predator-prey equations over a period of 5 years. The program should make a plot of the results.
Run your program for STEP_NUM = 100, and 1,000 and 10,000. We expect the correctness of the program to improve as we increase the number of steps, so believe your later results more than your first ones.
Turn in plots for STEP_NUM = 1,000 and for STEP_NUM = 10,000.
Question 3: HEAT, a boundary value problem.
In a boundary value problem, we are told the value of some quantity on the surface or edges of a region, and we are told the rules for how the quantity varies inside. This is sometimes enough information to determine the value of the quantity everywhere.
The simplest boundary value problem can be posed over a straight line, with known values at the two endpoints. We will look at a version of the heat equation. In this case, we are studying a metal rod, say. We also assume that parts of the rod are heated or cooled, perhaps by torches or fans. Whatever is going on, however, has been going on so long that the temperatures in the rod have stopped changing. They may vary from spot to spot, but they don't change over time. This is called the steady heat equation.
The main thing about heat is that it spreads out. More precisely, at any point X(I), the temperature T(i) "wants" to be equal to the average of the values at neighboring points. For physical reasons, the right quantity to look at is the negavite second derivative:
- d/dx d/dx T(x) = F(X)
If F(X) = 0 (no internal heat sources) then this says that the
temperature inside the rod is the straight line function between
the endpoint values.
To get a finite difference estimate of the second derivative, we essentially evaluate the first difference at T(I+dX/2) and T(I-dX/2):
d/dx d/dx T(x) = approximately
( T(I+1) - T(I) ) - ( T(I) - T(I-1) )
----------------- -----------------
DX DX
- -----------------------------------------
DX
which simplifies to
- T(I+1) / DX / DX
+ 2 * T(I) / DX / DX
- T(I-1) / DX / DX
and so our approximate steady heat equation can be written as
- T(I-1) + 2 * T(I) - T(I+1) = DX * DX * F(X(I))
Because each value T(I) is associated with two other unknown values, these equations are harder to solve than the initial value problems we saw before. Luckily, they are not very much harder. If we have set up N equally spaced points on the rod, including the beginning and end, then we have N unknowns and N equations (namely, the N-2 discrete heat equations, and the two boundary conditions.)
Let's write these equations as a linear system:
#1: T(1) = T1
#2 - T(1) + 2 * T(2) - T(3) = DX * DX * F(X(2))
#3: - T(2) + 2 * T(3) - T(4) = DX * DX * F(X(3))
#4: - T(3) + 2 * T(4) - T(5) = DX * DX * F(X(4))
... ... ...
#N-1: - T(N-2) + 2 * T(N-1) - T(N) = DX * DX * F(X(N-1))
#N: T(N) = TN
So all we need to do is store the coefficient matrix as A, the right hand side as rhs, and call MATLAB to solve this linear system.
Write a MATLAB program to solve the heat equation. Assume that the interval is 0 <= X <= 1.
A: solve the heat equation problem with a right hand side function F(X)=0. Assume that T(0.0) = T1 = 0.0 and that T(1.0) = TN = 100.0. Use N = 21 equally spaced nodes. Plot your solution.
B: solve the heat equation problem with a right hand side function F(X) which is 200 for 0.15 < X < 0.45, and 0 elsewhere. Assume that T(0.0) = T1 = T(1.0) = TN = 0.0. Use N = 21 equally spaced nodes. Plot your solution.