Dynamics

This chapter discusses the numerical computation of solu8tions of ordinary differential equations ma‘odeling dynamics systems. Differential equations are equations involving derivatives that wre used to determine a function, which is the solution. Here, we discuss ordinary differential equations, equations involving derivatives with respedt to a single variable, which we call t. In general discussions, we seek a function, x(t), and the differential equation is written in the form

ẋ = \frac{dx}{dt} = f(x) .  (1)

Here, both x and f have n components: x(t) = (x_{1}(t), ⋯, x_{n}(t)) ∈ R^{n}, and f(t) = (f_{1}(t), ⋯, f_{n}(t)) ∈ R^{n}. The components of x(t) determine the state of the system at time t, and f(x) represents the dynamics. In most real problems, f will depend on a number of problem parameters also.

For example, here is a model of a simple oscillator (such as a mass connected to a spring sliding on a table) subject to friction. The friction coefficient, γ also might change with time if friction causes the mass or the table to heat up. This is modeled through the equations

γ(t) = γ_{0} + γ_{1}θ(t) ,  (2)

where θ(t) is the temperature (say, in degrees centigrade) at time t and γ_{0} and γ_{1} are fixed coefficients. When γ_{1} > 0, the friction coefficient increases when the temperature increases. An oscillator with mass m, spring constant k, position q(t), and momentum p(t) is governed by the equations

mq̇ = p ,  (3)
ṗ = −kq − γq̇ ,  (4)

The −γq̇ term in the p equation represents the force of friction. The rate of work done by this friction force is force × q̇ = γq̇^{2}. We suppose the rate of heating is proportional to this rate of work and that heat dissipates with rate r, bringing the temperature back to the ambient temperature, θ_{0}. With this, the temperature evolution equation is

θ̇ = −r(θ − θ_{0}) + μγq̇^{2} .  (5)

We want to put this formulation into the generic form (1). The state of the system is determined by the three variables q, p, and θ. Although the friction coefficient, γ(t), is time dependent, it is determined by the other variables. For this reason, x(t) = (q(t), p(t), θ(t)). We identify f by writing explicit formulas for the components of x in terms of x and parameters. We find f_{1}(x) = ẋ_{1} by solving (3) for q̇:

q̇ = \frac{1}{m}p .  (6)

We find f_{2} = ṗ in terms of x, (i.e., in terms of q, p, and θ) by using (6) to put p instead of q̇ on the right, and by using (2) to eliminate γ in favor of components of x. This gives

ṗ = −kq − \frac{γ_{0} + γ_{1}θ}{m} p .  (7)

In the same way we can eliminate q̇ and γ from the right side of (??) to get

θ̇ = −r(theta − θ_{0}) + \frac{μ(γ_{0} − γ_{1}θ)}{m^{2}}p^{2} .  (8)

This puts the equations (3), (2), (4), (5) in the form (1), with

f(x) = ( \begin{matrix}\frac{1}{m}p \\ −kq − \frac{γ_{0}+γ_{1}θ}{m} p \\ −r(θ − θ_{0}) + \frac{μ(γ_{0}−γ_{1}θ)}{m^{2}}p^{2} \end{matrix} ) = ( \begin{matrix}\frac{1}{m}x_{2} \\ −kx_{1} − \frac{γ_{0}+γ_{1}θ}{m} x_{2} \\ −r(x_{3} − θ_{0}) + \frac{μ(γ_{0}−γ_{1}x_{3})}{m^{2}}x_{2}^{2} \end{matrix} ) .

Note that f depends on x and also on a large number of time independent parameters, m, k, γ_{0}, γ_{1}, r, and μ. This is typical.

We write the dynamics in the form (??) for the purpose of general mathematical discussion. In software it is fine to have something closer to the original formulation, such as:

int fEval( double x[3], double f[3]) { // Evaluate f for given x
. . . (parameter definitions)
double qDot, pDot, gamma, thetaDot;
double q, p, theta;
q     = x[0];  // Vector components start with 0 in C, as always.
p     = x[1];
theta = x[2];
qdot     = (1/m)*p;
gamma    = gamma0 + gamma1*theta;
pDot     = -k*q - gamma*qdot ;
thetaDot = -r*( theta - theta0 ) + mu*gamma*qDotqDot;
f[0] = qDot;
f[1] = pDot;
f[2] = thetaDot;
return 0; // Nothing can possibly go wrong.

The advantage of this point of view is that we can write a general program that solves any problem of the form (1).

We focus on the initial value problem. For this, in addition to the dynamics (1), we also specify initial conditions, x(t_{0}) = x_{0}. It is traditional to write x_{0} for the initial conditions without intending it to mean the zeroth component of x. In the initial value problem, x(t) represents the state of the system at time t. We presume to know the state at some specific time, t_{0}, and seek to compute the state at later times, t > t_{0}. The task, then, is to find x(t) for t > t_{0} given x(t_{0}).

Because of Newton’s law: F = ma, it is common in physical applications for problems to be formulated using second order differential equations. Suppose we have l particle coordinates q = (q_{1}, …, q_{l}) and corresponding force components F = (F_{1}(q), …, F_{l}(q)). Then Newton’s law is

m_{j}q̈_{j} = F_{j}(q) .

We put this into the standard format (1) using the extra variables p_{j} = m_{j}q_{j} (not the only possible choice). Then we have

q̇_{j} = \frac{1}{m_{j}} p_{j} ,
ṗ_{j} = F_{j}(q) .

This gives n = 2l, x = (q_{1}, …, q_{l}, p_{1}, …, p_{l}), and f_{j} = q̇_{j} = \frac{1}{m_{j}} p_{j}, and f_{j} + l = ṗ_{j} = F_{j}(q), for j = 1, …, l. In this process, we have replaced a system of l second order equations with an equivalent system of n = 2l first order equations.

We could reduce higher order differential equations as well. Consider, for example, the single third order equation \frac{d^{3}q}{dt^{3}} = q^{2}. For this, take x_{1} = q, x_{2} = q̇, and x_{2} = q̈. Then we have f_{1} = x_{2}, f_{2} = x_{3}, and f_{3} = x_{1}^{2}. A single third order equation has become a system of three first order equations.