These overheads taken from Informal DEs
"An equation between specified derivatives of an unknown function, its values, and known quantities and functions". For instance,
F(x,y,y') = 0
G(x,y,y',y'') = 0
a(x)y'+b(x)y + c(x) = 0
my'' + ky' +sy = g(x)
Our job is to find the function(s) y(x) so that our given equation holds.
Sometimes solution is possible through more or less difficult analytic techniques, sometimes we must resort to approximations, (numerical methods, numerical solutions), with which we're most concerned here.
A Vital Issue! Using physical laws to model physical systems as equations. Like algebra word problem. Setting up the equations is often 9/10 of the battle and we barely touch on it here.
Informally, the derivative y'(x) gives the slope of function y(x). Likewise the integral (say from 0 to t) of y(x) measures the area under the curve y(x).
The picture says it all:
The integral is the limit of a sum; it's a continuous sum of infinitesimals, same as a derivative is an infinitesimal difference.
In computing, we don't have infinitesimals nor continuous numbers: every "integral" is really a summation (unless it's a symbolic integral, maybe in Mathematica).
If we know y'(x) = g(x), then if we integrate both sides, we get an equation y(x) = G(x), where G is the integral of g.
Hence "Solving a DE" = "Integrating a DE" and "Numerical Integration" is approximate computer solution.
Examples:
d(ax^{n})/dx = nax^{n-1}, so
y' = nax^{n-1} means y = ax^{n} +C.
Here, C is a constant of integration. Generally DE solutions are families of equations.
sin(x)' = cos(x), cos(x)' = -sin(x)
d e^{ax}/dx = ae^{ax}, so y' = ay means, using t
instead of x this time,
y(t) = Ce^{at}.
As written a DE has a family of related solutions. For example,
If y'(x) = 20x^{3} + 14 x, then
y(x) = 5x^{4} + 7 x^{2} + C, where
C is a "constant of integration".
Can specify one of a family of solutions by specifying Boundary
Conditions (values of variables and their derivatives) it must also satisfy.
E.g. for a DE like
y'(t) = 3t^{2} + 13, the boundary condition
y(0) = 19 identifies a unique solution:
integrating, get y(t) = t^{3} +13t + C, so must have
C = 19 to get y(0) = 19.
We'll only have Initial Conditions (hence "Initial Value Problems"). Another case we won't explore is the one where we need to specify what happens at "both ends" of the solution for instance ("Boundary Value Problems"). E.g. What's the distribution of temperature in a metal plate whose periphery is kept cooled to a constant temp (the boundary condition) but which is being heated in the middle?
Ordinary and Partial DEs: ordinary d/dt and partial ∂ / ∂ t derivatives.
Order of DE: order of highest derivative appearing in it. An nth order DE can be expressed as a system of n 1st-order DEs, as we'll see.
Linear and Nonlinear DEs: Linear DE in y is a sum of weighted
y-derivatives of different orders set equal to 0:
a(x)y'' +b(x)y' + c(x) y + d(x) = 0
Nonlinear has more complex relations between the derivatives, like
a(x)y''y +b(x)(y'^{2}) + c(x) √ y + d(x) = 0
Linear constant coefficient equations have solutions in exponentials
(decays, explosions, and sines)
ay'' +by' + cy + d = 0
Homogeneous and Nonhomogeneous: Hom looks like
a(x)y'' +b(x)y' + c(x) y = 0
Nonhom:
a(x)y'' +b(x)y' + c(x) y = -d(x)
Homogeneous can be a steady state solution (no 'driving force' from outside). Any solution to a DE can have the homogeneous soln added to it and stay a soln.
Cute trick: change a higher-order ODE (higher-order derivatives) to a set of related first-order ODEs (first derivatives only).
Invent new variables (say u_{1}, u_{2}, u_{3},
etc.) so that the first one is y', the derivative of y; the second is the
derivative of the first, (y''); etc. So if
y''' + g(x)y'' +r(x)y' + s(x)y = t(x),
then we must have
u_{1}(x) = y'
u_{2}(x) = u_{1}' = y''
u_{3}(x)= u_{2}' =
y''' =
t(x) -s(x)y -r(x)u_{1} -
g(x)u_{2}.
There is a rather nicer (to my mind) matrix form of this, also given in the pdf reading: Informal DE Intro.
y''' + g(x)y'' + r(x) y' + s(x)y = t(x)
|y | | 0 1 0 ||y | |0 | d | | | || | | | -- |y' |=| 0 0 1 ||y' |+|0 | dx | | | || | | | |y''| |-s(x) -r(x) -g(x)||y''| |t(x)|
Another example of a nonlinear system of three first order DEs (The Lorenz
Attractor):
x' = 10(x-y)
y' = -xz +28x -y
z' = xy -8z/3
Which has a "chaotic" regime, solutions as below: red and blue solns from almost identical initial conditions.
Final word. All the equations we feed Matlab to get our ODE
solving done will look like
y' = f(y,x).
For Initial Boundary Value Problem. Assume know y(0), the given initial condition. But y'(x) tells you how much to change y for a small step in x. So step a bit dx along x and use y' to get a new value for y: in fact, just iteratively set y(x+dx) = y(x) + (dx)y'.
That's ALL THERE IS!! Everything else is just bells and whistles.
Add Rectangles, Trapezoids, or small areas from higher-order approximations: Trapezoidal Rule, Simpson's Rule, Runga-Kutta, Predictor-Corrector... All doing the same idea with varying degrees of accuracy, speed, flexibility, efficiency.
Stiffness: solutions at hugely different timescales means intelligent and varying choice of step size is needed (e.g. Matlab's ode45 and ode45s).
Insult Warning! This stuff is simple!
Here's an example (PLENTY more out there in Matlab documentation , etc.
Lorenz Attractor:
w' = -xz +28x -w
z' = xw -8z/3
First Matlab needs a function that basically copies the equations above to
construct a column vector
of derivatives
y' = f(t,y)
with one row
for each of the unknowns y(i)
in our equation.
In our equations there is only dependence on x,w,z, not on t.
Our derivatives
function takes two
arguments, a scalar 'time' t
and a vector of unknown
variables.
So we need a little renaming of variables in the equations:
y(1) = x
y(2) = w
y(3) = z
Then we write directly from the equations (using the new
y()
names) a function that looks like this:
function dy = lorenz(t, y)
% Lorenz Attractor. t a scalar,
% y a vector of y(i)'s.
% dy returned as a column vector
% of dy(i)'s
% the orig x,w,z are
% now y(1), y(2), y(3)
dy = zeros(3,1);
dy(1) = 10*(y(2) - y(1));
dy(2) = - y(1)*y(3) + 28*y(1) - y(2);
dy(3) = y(1)*y(2)-(8/3)*y(3);
end
The equations are implicitly parameterized by time, which steps along according to a timespan parameter vector (next slide).
Now we've got our derivatives function, which computes the current derivative value from the values of lower-order derivatives and the time.
We want to pass this function to another function as argument. As usual
with imperative languages, we'll pass a pointer, or a
handle for this function, given by Matlab's @
operator: E.g., either the bare
@lorenz
Or we can give it a name:
lorenzhandle = @lorenz
We need to
specify the time interval to be
covered by the integration. The documentation says it can be a 2-vector
[start, finish]
or an N-vector like
timespan = linspace(start, finish,1000)
With the 2-vector, Matlab makes up the spacing: the N-vector
provides explicit times (not necessarily evenly spaced).
We also need a (row or column) vector of initial conditions for the
system.
E.g. for Lorenz,
let's start at some point in space
Let's run for 20 time-units divided up into 1500 steps. We
also must be able to spell ode23, the simplest and cheapest
of Matlab's many ODE solvers.
We expect the output
to be a column vector of times (the same as our timespan input, or
made up by Matlab to lie between our start and end time), and
values of the function (and any lower-order derivatives it computes
while integrating). OK, we're off!...
%Script to Run Lorenz
lorenzhandle = @lorenz;
y0 = [10, -10, 20]; % initial conditions
%next specify 1500 instants
timespan = linspace(0, 20, 1500);
% ** boom **
[T,Y] = ode23 (lorenzhandle,timespan, y0);
% ** it's over! **
plot(Y(:,1),Y(:,3), 'r');
hold on
% new initial conditions close to previous
y0 = [10, -10, 20.1];
[T,Y] = ode23 (lorenzhandle,timespan, y0);
%BUT big difference in the two solns
plot(Y(:,1),Y(:,3), 'b');
In a vacuum, a falling body accelerates due the force of gravity.
On earth, the acceleration due to gravity is called
Let's plot the velocity and position of a ball (in a vacuum) acted on by gravity. We can start the ball at what ever height and velocity we want... throwing it up sounds interesting (and is like the ballistics assignment).
The equation governing y as a function of time says that the
acceleration of a body due to gravity is in the "down" direction and
is always equal to g.
Assuming
y'' = -g
That's it...that's our DE to solve. Notice it's 2nd order so we'll
make a system of 2 simpler ODES.
We put that into our favorite matrix form this way.
d |y | |0 1| |y | |0 | -- | | = | | | | + | | dx |y'| |0 0| |y' | |-9.8|
or:
dy/dx = y'
d^{2} y / dx^{2} = y'' = -9.8
So our derivatives function in Matlab is this:
function dy = ball_der(t,y)
dy = zeros(2,1);
dy(1) = y(2); %dy/dx = y'
dy(2) = -9.8; %d^{2}y/dx^{2}=y''=-9.8
end
We now get to choose a start and finish time and initial conditions. May as well start at zero time and altitude. Let's start the ball off with an upward velocity of 100m/sec.
% script
ballhandle = @ball_der;
tspan = [0 20]; % [st fin] (secs)
y0 = [0 100]; % [ y y'] initial
[T Y] = ode23(ballhandle, tspan, y0); %boom
....
>> [T Y]
>> plot(T,Y(:,1))
>> plot(T,Y(:,2))
Notice Matlab starts out with small
T Y VEL 0 0 100.0000 0.0000 0.0001 100.0000 0.0000 0.0005 100.0000 0.0000 0.0025 99.9998 0.0001 0.0125 99.9988 0.0006 0.0625 99.9939 0.0031 0.3124 99.9694 0.0156 1.5613 99.8469 0.0781 7.7826 99.2344 0.3906 38.3148 96.1719 1.9531 176.6205 80.8594 3.9531 318.7392 61.2594 5.9531 421.6580 41.6594 7.9531 485.3767 22.0594 9.9531 509.8955 2.4594 11.9531 495.2142 -17.1406 13.9531 441.3330 -36.7406 15.9531 348.2517 -56.3406 17.9531 215.9705 -75.9406 20.0000 40.0000 -96.0000
Velocity:
Height:
Simulations are great fun, since we get to play around with the rules. Suppose the force of gravity were a nice even 10, and that its direction is reversed 10 seconds into the ball's flight?
function dy = gball_der(t,y)
if t<10
g = -10;
else
g = 10;
end
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = g;
end
and
Height:
Our goals were:
Trajectory of human head launched at 40 degrees of elevation into enemy redoubt. Flight took about 15 seconds. Distance unit is the meter (about 3 feet for you Angles and Saxons).