Copyright © Michael Richmond.
This work is licensed under a Creative Commons License.
See Sections 1.1, 25.1-25.2 of your textbook for more information on these methods.
A differential equation is one which expresses the change in one quantity in terms of others. Thus,
2
y = x
is not a diff eqn, whereas
dy
-- = 2x
dx
is a diff eqn.
Diff eqns occur very frequently in all branches of physics, and so we must devise ways to deal with them. Here are several:
Sometimes, when you're lucky, the physical situation has an analytic solution. For example, consider the relationship between velocity and time, for the case of constant acceleration:
dv
-- = a = -g
dt
v = at + v0 = -gt + v0
Pretty simple! In this case, we can even find an exact relationship
between velocity, position and time:
dy
-- = v = -gt
dt
1 2
y = - --- g t + y0
2
If you can find an exact analytic solution, use it. That's always the best way to go, since you can figure out the quantities of interest (position, velocity) for any time, with absolutely no error.
Unfortunately, in many situations, the physical relationships between quantities is more complicated, such as
dv 2
-- = a = -g + K v
dt
In most of these cases, there is no exact solution --
or, if there is one, it might be very hard to figure out.
Under these circumstances, we can apply
numerical methods to understand the behavior of the
system -- but we always fall back on analytic solutions
to check the numerical results.
Let's consider the following real-life problem: Joe Geronimo jumps out of an airplane at an altitude of H = 100 meters. What are his altitude y(t) and velocity v(t) as functions of time? How long until he hits the ground?
We'll start by making a simplification: assume that there is no air resistance. We'll include air resistance and other complications in the future.
In this simplified situation,
2
d y dv
------ = ---- = a = -g = -9.8 m/s^2
2 dt
dt
This can be solved analytically!
v(t) = v0 - g*t
where v0 is a constant of integration.
We use the initial conditions
of the problem to set its value.
In this case, Joe's speed at time t = 0, when he
jumps out of the plane, is v = 0 m/s.
So that means v0 = 0 m/s, and so
v(t) = 0 - g*t
Likewise, we can solve for Joe's position as a function of time:
y(t) = y0 + v0*t - 0.5*g*t^2
Once again, there are constants of integration.
We know that at the start, Joe's height is H = 100 m,
and his velocity is v = 0 m/s, so we have
y(t) = H + 0 - 0.5*g*t^2
Okay, that's the analytic solution. We could use it to solve for Joe's speed and altitude at any time, and we could also figure out the time at which Joe's altitude is zero (ouch!).
Now, let's see how we can solve the problem numerically. Consider Joe's velocity first. We know
dv
-- = -g
dt
and we want to find v(t).
Hmmm. Suppose that we know Joe's velocity at some particular instant t0. As we can for any arbitrary function, we can expand his velocity as a function of time in a Taylor series:
1 2
v(t) = v(t0) + v'*(t - t0) + - v''*(t - t0) + ....
2
In this particular case, we happen to know the first derivative of velocity, v'. If we ignore all the terms of second, third, and higher order, we can calculate Joe's speed an instant later:
v(t) = v(t0) + v'*(t - t0)
Note that this solution only gives the correct answer if the second, third and higher-order derivatives are all zero. However, if we pick a small enough time interval, then it will give us a result which is close to the true answer, even without the terms of higher order.
When dealing with numerical solutions, we often write quantities with little subscripts to indicate "current" and "next" values:
t = current time
i
t = next time
i+1
v = current velocity
i
v = velocity at next time
i+1
v' = current derivative of velocity
i
In which case, we can write:
v = v + v' * (t - t )
i+1 i i i+1 i
In many cases, it helps to define a "timestep:" the amount of time which passes between each of the moments at which we calculate position, velocity, or other physical quantities.
dt = t - t
i+1
Code to solve differential equations almost always contains a loop somewhere in which quantities at the current time are used to calculate quantities at the next instant:
LOOP
know t v v'
i i i
v = v + v' * dt
i+1 i i
v' = (something)
i+1
t = t + dt
i+1 i
% now prepare to go through loop for the next moment
v = v
i i+1
v' = v'
i i+1
t = t
i i+1
END_OF_LOOP
This approach is called Euler's method: it uses a first-order approximation to the Taylor series to calculate new values for physical quantities.
Now, in this particular case, as Joe falls without air resistance, we know
dv
v' = -- = -g
dt
2
d v
v'' = ---- = 0
2
dt
and all higher derivatives of velocity are also zero.
So the Taylor series turns into
v(t) = v(t0) + (-g)*(t - t0) + 0 + 0 ...
v(t) = v(t0) + (-g)*(t - t0)
which means that Euler's method
v = v + (-g) * (t - t )
i+1 i i+1 i
will give the exact value for velocity.
Lucky us!
Most of the time, the second (and higher) order term of
the Taylor series is not zero, and so Euler's method
yields only approximately correct values.
For example, consider the height of Joe above the ground as a function of time: y(t). We can write a Taylor series for it:
1 2
y(t) = y(t0) + y'*(t - t0) + - y''*(t - t0) + ....
2
Now, in this case,
y' = v
y'' = a = -g
y''' = 0 (and all higher derivs, too)
The second derivative of position with respect to time is not zero (in fact, it's a constant, -9.8 meters per second squared). So, if we apply Euler's method to calculate position as a function of time:
y = y + v * (t - t )
i+1 i i+1 i
Let's try that and see what happens:
t a v y (s) (m/s^2) (m/s) (m) ----------------------------------------------- 0 -9.8 0 100 1 -9.8 -9.8 2 -9.8 -19.6 3 -9.8 -29.4 4 -9.8 -39.2 -----------------------------------------------
Q: If we use this scheme,
how far above the ground is Joe at t = 4 s?
Q: Estimate when Joe hits the ground.
But this is not the right answer: the analytic solution indicates that Joe reaches the ground at t = 4.52 s. Uh-oh! What can we do!? Well, we have several options:
1 2
error = - * g * (t - t )
2 i+1 i
If we make the timestep small, then the error --
which goes as the square of the timestep -- will become
very small; small enough, perhaps, that it will make
no significant difference in our result.
1 2
y = y + v * (t - t ) - --- * g * (t - t )
i+1 i i+1 i 2 i+1 i
In this case, because the second derivative is constant
and all higher derivatives are zero, this ought to give us the exact
value for the next position.
In general, however, the second- and higher-order derivatives
are not so simple. Often, the higher derivatives change
with time.
This week, let's stick with Euler's method, for both velocity and position, to solve the problem. The initial conditions are:
dv
---- = a = -g = -9.8 m/s^2
dt
Copyright © Michael Richmond.
This work is licensed under a Creative Commons License.