Differential Equations
Contents
Differentiation
Integration
Linear Differential Equations
Nonlinear Differential Equations
Numerical Solution for 2nd Order Differential Equations
Contour Integrals
Derivation of some Differential Formulas
Following is a very simple differential equation applicable to exponential decay such as the decay of radioactive element:
dy(t) / dt = -k y(t) ---------- (1)
where the quantity of the material y(t) is a function of the time t only, dy(t) = y(t2) - y(t1) is an infinitesimal change of y in an infinitesimal time interval dt = t2 - t1 with t2 > t1. In plain language, it states that the reduction of y, i.e., "dy" (and hence the minus sign) over a time interval "dt" depends on the remaining amount of y. The parameter k is a proportional constant which is called the rate of decay in this case, and is unique to a particular system. In particular, if k = 0, then y is not decaying at all.
The concept of differentiation depends critically on the fact that a small number such as 0.00000006 divided by another small number such as 0.00000003 yields a finite value 2 (in this example). Thus dy and dt may each be infinitesimal, but the derivative dy/dt is finite. Graphically, dy/dt measures the slope of a curve as shown in Figure 02. The "d" in calculus always signifies an infinitesimal change.
[Top]
ln(y) - ln(A) = -kt, or y = A e-kt
where e = (1 + 1/n)n 
2.71828... as n 
is the exponential base used to simplify the form of many formulas (it plays a special role in calculus similar to the "
" in geometry), "ln" is the logarithms in this base, and A is the integration constant (which vanishes and the original differential equation is recovered if we reverse the process by taking the differentiation again) to be determined by the initial condition. For example, if at t = 0, y = yo, then A = yo. Alternatively, the integration constant can be replaced by the definite integral such that 
f(y) dy = F(y)
= F(b) - F(a), where F(y) is the result of the integration.
 |
In either way, the final solution for Eq.(2) is:
y = yo e-kt --------------------------------- (4)
The green curve in Figure 02 shows the radioactive decay of Iodine 131 with a half life of 8.1 days. Half life is the interval when half of the original is gone. This graph corresponds to a value of 1/k = 8.1/loge(2) = 11.57 days, and yo = 1 gm. in Eq.(4).
|
|
|
Alternatively, Eq.(3) can be solved numerically by a computer. If we denote 
y to be the width of the strips in Figure 01, then the area of the strip y/y = -kt. The relationship between t and y can be approximated by the summation (together with the initial condition at t = 0, y = yo):
ny/y = -k nt = -k [(0 - t1) + (t1 - t2) + 
+ (tn-1 - tn)] = k tn ---------- (5)
where n denotes the number of strips, and the vaule of y can be taken at the middle of each one (as y1, y2, , yn). The result (for a given y) becomes progressively more accurate as the width of the strip is reduced at the expense of computing time. Numerical integration is the only option to solve the differential equation when there is no analytical solution. Eq.(5) is a very simple example, numerial solution can be very complicated involving iterations and finding the root of an equation.
[Top]
In the more general case when the function depends on more than one variable, the differentiation would have to be performed against one variable at a time while keeping the others fixed. For example, the partial derivative of x for f(x,y) is defined by:
 | ---------- (6) |
Many differential equations in physics involve multiple variables. Such kind of equations are referred as partial differential equations. The wave equation in three dimensional space is an example:
 | ---------- (7) |
where x, y, z are the spatial coordinates, t is the time, v is the propagation velocity, u(x,y,z,t) is the displacement of the wave, and the second derivative is defined as d2u/dx2 = d(du/dx)/dx. For wave motion in one dimension, e.g., along the x axis, the 2nd and 3rd terms in Eq.(7) vanish; the solution for u at any point x and any time t is given by the cosine (or sine) function:
 |
u(x) = A cos(kx -  t) ---------- (8)
where = 2 , k = 2/ , and A is the amplitude. The frequency is related to the wavelength and wave velocity by the formula: = v. Figure 03 shows the cosine wave (for A = 1) at two instants of t = 0, and t = /.
|
|
Note: see footnote  for derviation of the differentiation for the sine and cosine functions.
|
Eqs.(1) and (7) are linear differential equations because they involve only linear terms in either y(t) or u(x,y,z,t).
Examples can be found in many physical systems such as the harmonic motion in "Classical Mechanics", and the Schrodinger equation in "Quantum Mechanics". Linear equation (such as Eq.(1) or Eq.(7)) has the special property that if y1 and y2 are the solutions of k equals to k1 and k2 respectively, the combined solution y = A1y1 + A2y2 is again a solution of the same differential equation, where A1 and A2 are arbitrary constants. This is known as superposition in wave theory such as the one described by Eq.(7). Superposition forms the base in explaining interference and the production of wave packets. The mathematical formula for a wave packet is just the generalization of the superposition with the wave solution in the form of Eq.(8) (at t = 0 for example):
 | ---------- (9) |
 |
where the coefficients An's are replaced by the weighting function f(k - ko), and the sum becomes integration over a continuous range of k from - to + (see Figure 04a). At x = 0 the cosine function is equal to 1, hence all contributions to the integral comes from the weighting function, and the result is large. As x increases, the cosine function becomes a rapidly oscillating function of k, and its integral tends to cancel out resulting in the form of a pulse (wave packet) as shown in Figure 04a. In the limiting case where the weighting function is expressed in term of a delta function - A (k - ko), the wave packet is reduced to the monochromatic form such as shown in Figaure 03 or Eq.(8) with k = ko; the wave is now extending from - to +. In general, there is a certain spreads x in x, and k in k. In quantum theory, the uncertainty principle postulates that xp >  =1.054x10-27erg-sec,
|
|
where the momentum p = k comes from the de Broglie formula linking p and . Note:
It can be shown that the inverse of Eq.(9) is given by f(k - ko) = (1/) y(x) coskx dx, |
which together with Eq.(9) are known as Fourier transforms. Similar linear combination with the exponential function in Eq.(4) is called Laplace transforms. The Laplace transformation offers a direct method of solving linear ordinary and partial differential equations (the mathematical detail will not be pursued here).
By virtue of the orthogonality relation:
|
 , |
 |
the cosine function in Eq.(9) acts like an unit vector1 in the 3-dimensional space (Figure 05), such that inim = nm where i1 = i, i2 = j, and i3 = k, and nm = 1 for n = m; = 0 otherwise. Thus, in quantum mechanics we can consider the wave function in force free space (like the wave packet) as an expansion in an infinite dimensional space with cos(nx) as its bases (unit vectors). Similar expansion can be expressed with wave function in bound states (where n and m become integers) for the cases of non-vanishing potential V(x). The idea has been expanded even further in quantum field theory to the more abstract Fock or Hilbert space, where it is just a mathematical symbol such as |nk , where nk represents the number of particles in the state k with certain energy-momentum, spin, ... etc. The
|
Figure 05 Unit Vectors (The vector v can be expressed in terms of the unit vectors)
[view large image] |
corresponding orthogonality relation is now expressed as  nk|nk' = nknk'. In quantum theory, the coefficient or component (such as the f(k-k0) in Eq.(9)) for each of such basis has been identified as the probability amplitude with the state k.
|
1A vector is a mathematical entity with two or more space dimensions such as the velocity of an object v in Figure 05. Loosely speaking, a (second rank) tensor is the product of two vector spaces (un-prime and prime) with ii', ij', ..., kk' as its bases.
[Top]
In a nonlinear differential equation such as:
dy(t) / dt = -k' y2(t) ---------- (10)
with solution in the form:
1/y - 1/yo = k't or y = 1 / [k't+(1/yo)] ---------- (11)
(see formula for differentiation of function with power n in footnote¶), it is obvious that the combined solution y = A1y1 + A2y2 is not a solution of the same differential equation any more. The blue line in Figure 02 is the decay curve described by Eq.(10) with k' = -(dy/dt)/y2|t=0 ~ 0.07/day-gm. The red line represents the case with a constant decay rate, i.e., dy/dt = -ko. Comparison shows that the negative feedback from y increasingly slows down the process as the contribution changed from y0, to y1, and y2.
 |
The logistic equation for the growth of population is a little bit more complicated but is still solvable in closed form:
dy/dt = ry(1 - y/K) ---------- (12)
where y is the population size, r refers to the rate of reproduction, and K denotes the carrying capacity (the maximum size of population allowable by the environment). Following similar procedure for solving Eq.(2) in footnote §, we can derive the integration in the form:
 dy/y(1-y/K) = d{ln[y/(1-y/K)]} = rdt ---------- (13) |
|
which yields: y/(1-y/K) = [yo/(1-yo/K)] ert or y = Kyoert/[K + yo(ert - 1)] ---------- (14)
|
In the limit t , the population size converges to the equilibrium x = K; while for the case K (unlimited population size), it reverts back to exponential growth as opposed to the exponential decay in Eq.(2) (see Figure 06). This is an example to demonstrate clearly the effect of negative feedback in controlling runaway process.
If we redefine the meaning of the parameter such that r = (a - b) and re-label the variable y with y + x = 1, Eq.(12) can be used to describe the process of natural selection, where a and b are the reproductive rates for species A and B, y and x are the corresponding population numbers for K = 1. Without further ado, it is clear that if a > b, then A will become more abundant than B. Eventually A will take over the entire population; B will become extinct (and vice versa).
 |
Eq.(12) can also be recast to the form of difference equation with discrete steps, e.g., using the generation number n as running variable with n = 1 instead of the time t with t0:
yn+1 = Ryn (1 - yn) ---------- (15)
where R = 1 + r, and K = 1, which rescales the population size y to between 0 and 1. Thus, |
|
R < 1 means r < 0, which implies decay instead of growth. The behavior of Eq.(15) is more complicated depending on the value of R as shown in Figure 07 and summarized below (for yo = 0.1, and 80 generations): |
- For 0 < R < 1 - The population will go extinct.
- For 1 < R < 3 - The equilibrium population size is given by y* = (R-1)/R (by setting yn+1 = yn in Eq.(15)). All trajectories starting from any initial yo will converge to this value. The point y* is called global attractor.
- For 3 < R < 4 - The equilibrium y* becomes unstable. For values of R slightly above 3, we find a stable oscillation of period 2. As R increases, the period 2 oscillator is replaced by period 4, then by 8, and so on. For R = 3.57 there are infinitely many even periods. For R = 3.6786 the first odd periods appear. For 3.82 < R < 4 all periods occur.
- For R = 4 - The dynamics becomes unpredictable in the sense that the trajectories starting from two slightly different initial points will diverge completely. This is called deterministic chaos, and is one of the characteristic of nonlinear equation(s).
- For R > 4 - Values of y > 1 can be generated, and therefore negative y in the following year, thus the equation becomes unphysical (unbiological) beyond this point.
As shown in Figures 05 and 06, the behavior of the logistic differential and difference equations is drastically at odds with each other. The crucial difference is the absence of the k parameter, which can keep minimize the negative feed back. Nonlinear differential equations appear in "General Relativity", "Fluid Dynamics", ..., and the Lorenz systems (in chaos theory), which is described by a set of coupled differential equations:
dx/dt = 
(y - x), dy/dt = rx - y -xz, dz/dt = xy - bz ---------- (16)
where t is time (the independent variable), x, y, z are the dependent variables, and , r, b are the parameters. A characteristic of nonlinear differential equation(s) is the mixing of many functions (dependent variables) or the function acting on itself. While it seems to be more realistic, the solution for such holistic approach is more difficult to find. Sometimes the nonlinear equations can be linearized because some variables (or their derivatives) are small in comparison with the others. See Chaos Theory for more about nonlinear equations.
[Top]
In general, a 2nd order differential equation involves the independent variable x, the dependent variable y, its 1st derivative dy/dx, and the 2nd derivative d(dy/dx)/dx. Numerical solution requires only simple programming whenever we can isolate the 2nd derivative in the form:
d(dy/dx)/dx=f(y',y,x),
where y'=dy/dx=z.
Thus, we can rewrite the original equation as:
dz/dx=f(z,y,x).
In the form of difference equation:
z(n+1)=z(n)+f[z(n),y(n),x(n)]dx
y(n+1)=y(n)+z(n)dx
Input parameters include the step size dx, the number of steps max, and and the initial conditions x(0), y(0), z(0).
Numerical integration is carried out by iteration as shown below (in Basic, also see actual example for a cosmological model):
input max
input dx
input x0
input y0
input z0
n=0
zn=z0
yn=y0
xn=x0
#from the actual expression of the function at x0
fn=f0
MainLoop:
znn=zn+fn*dx
ynn=yn+zn*dx
xnn=xn+dx
n=n+1
if (n-max)=0 then goto Ende
zn=znn
yn=ynn
xn=xnn
#from the actual expression of the function
fn=fnn
print xnn;:print ", ";:print ynn
refresh
goto MainLoop
Ende:
The numerical solution for 1st order differential equation follows the same idea only simpler (see a Basic Program for the Logistic Equation). Other "home-made" programs such as Base Conversion (between decimal, hexadecimal, and binary), Unit Conversion (fps to SI, SI to fps), Mathematical Functions, and Loan Payment (all in Basic) are available (just clicking the link) to show the subtlety of computer programming.
[Top]
 |
It was known in the 17th century that the roots of polynomial equations, i.e., when f(x) = 0 in:
f(x) = anxn + an-1xn-1 + +a1x + a0 (such as f(x) = x2 + 2 = 0 at x =   ),
often involve square roots of negative numbers. It was realized subsequently that the complex number, which is defined as the combination of a real number and a square root of negative number, is more useful than the real number alone. In the Cartesian form the complex number z is expressed as z = x + iy, where the negative sign in the square root has been absorbed into the symbol i =  . In terms of the polar coordinates r and  , z = r ei, where r2 = x2 + y2, and tan = y/x (Figure 08). The complex conjugate of z is defined as z* = x - iy or z* = r e-i, such that z z* = r2. The term "imaginary" for the part associated with "i" was coined by |
|
René Descartes and was meant to be derogatory. The study of functions of a complex variable is known as complex analysis. Unlike real functions, which are commonly represented as two-dimensional graphs, complex functions have four-dimensional graphs.
|
 |
Contour integral is related to a special class of complex function called analytic function denoted here by f(z) = u + iv, with the complex variable z = x + iy. It has the property that the derivative df(z)/dz exists uniquely at every point. In other words, the way z 0 is indepen-dent of the path, e.g., we may choose a path with x=0 or y=0. It follows that the real and imaginary parts of an analytic function are solutions of the two-dimensional Laplace equation: |
|
 ,  .
|
It can be shown that the contour integral of an analytic function f(z) around any closed path c is zero, i.e.,
 | ---------- (18) |
This is illustrated in diagram (a) of Figure 09. The vanishing integral is related to the fact that the number of lines entering and leaving the contour are the same. In case there is a source or sink (collectively called a pole) within the region so that the number of lines entering and leaving the contour is not balanced as illustrated in diagram (b) of Figure 09, the contour integral becomes:
 | ---------- (20) |
where z0 denotes the position of the pole. For example,
 | ---------- (21) |
In general, for a function with multi-poles such as f(z) / (z-z1)(z-z2)...(z-zn), where f(z) is an analytic function, the contour integral enclosing all these poles is given by 2i [f(z1)/(z1-z2)...(z1-zn) + f(z2)/(z2-z1)...(z2-zn) + ... + f(zn)/(zn-z1)(zn-z2)...].
A technique to evaluate ordinary integration is to use the contour integral with poles whenever Eq.(17) is satisfied. For example, considering the following integral:
 |
I =  dx / (x2 + 1) ---------- (22).
This integral is part of a closed contour integral, which has a pole at z = i (see Figure 10):
Ic = dx / (x2 + 1) + cR dz / (z2 + 1) = I = 2i [1/(z + i)]z=i = 2i (1/2i) = --- (23), |
|
where the second term vanishes if we let the radius of the semi-circle R .
|
Another neat trick is to create pole(s) inside the loop (contour). Considering the following integral:
I = dx / (1 - x2) ---------- (24)
There are two poles at x = +1 and x = -1, right on the real axis. In order to apply the loop in Figure 10 to evaluate this integral we can artificially move one of the poles a wee bit up and the other one down by an amount i
as shown in Figure 10, and then let 0 at the end. Thus,
Ic = dx / (1 - x2) + cR dz / (1 - z2) = I = 2i {[1/(1+ i+z)]z=1+i}0 = 2i (1/2) = i ---------- (25).
Actually, the answer can be obtained as easily by substituting y2 = -x2, dx = i dy in Eq.(24), and using the result from Eq.(23).
[Top]
¶ Differentiation of power function (for n 
0, see Eq.(32) for the case with n = 0):
d(yn) 
(y + y)n - yn yn (1 + y/y)n - yn 
yn (1 + ny/y) - yn nyn-1y = nyn-1dy ---------- (26)
where we have used the power series expansion: (1 + x)n = 1 + nx + n(n-1)x2/2! + . In Eq.(10), we can re-arrange the y variable to one side and re-write y-2dy = -d(y-1) by applying the formula in Eq.(26) for the case of n = -1. Integration to obtain the solution becomes trivial once it is written in this form.
Note: Only the highest term or the 2nd highest term (in case of cancellation of the highest term) is taken in the series expansion. Readers of this page are strongly urged to complete the intermediate steps for many of the derivations as it offers a chance to apply their skill on the subjects.
Derviation of the differentiation for the sine and cosine functions:
dsinx = sin(x + x) - sinx = sinx cosx + cosx sinx - sinx 
cosx x = cosx dx ---------- (27)
dcosx = cos(x + x) - cosx = cosx cosx - sinx sinx - cosx -sinx x = -sinx dx ---------- (28)
where we have used the series expansion (derived from the Taylor series expansion about the point xo = 0):
sinx = x - x3/3! + ---------- (29a)
cosx = 1 - x2/2! + ----------(29b)
and the trigonometric formulas:
sin(A + B) = sinAcosB + cosAsinB ---------- (29c)
cos(A + B) = cosAcosB - sinAsinB ---------- (29d)
These trigonometric formulas were discovered more than one thousand years ago by the Muslim mathematician Abu Wafa.
§ Derivation of the relation d[ln(y)] = dy/y:
 |
Let us start from the differentiation of the logarithmic function loga(y) with an arbitrary base a:
d[loga(y)] loga(y + y) - loga(y) loga(1 + y/y)
(y/y) loga(1 + y/y)(y/y) ---------- (30)
where we have applied the logarithmic relations (Figure 11):
loga(u/v) = loga(u) - loga(v) ---------- (31a)
loga(u)n = n loga(u) ---------- (31b).
|
|
|
Now if we define n = y/y, e
(1 + 1/n)n, a = e, ln(y) = loge(y), and use the identity loga(a) = 1; then
d[ln(y)] y/y = dy/y -------- (32).
As a corollary from the definition of the exponential base e, it can be shown that differentiation of the exponential function ex is simply:
d(ex) = ex dx ---------- (33).
Incidentally, logarithm is not equal to "log-a-rhythm" (Figure 12). It was invented for manipulating large number by reducing it to a manageable size, e.g., log10(1000000000) = log10(109) = 9. The first logarithmic table was constructed by the Scottish mathematician John Napier (1550-1617) in the base of e (since known as natural log). Also note that in term of complex function with a real and an imaginary parts, the exponential and trigonometric functions are related by the formula:
 |
y = ei x = cos x + i sin x ---------- (34)
This relationship was discovered by the Swiss mathematician Leonhard Euler (1707-83). The modern derivation can be obtained easily if we perform a Taylor series expansion about the point xo = 0 for ex, substituting to the independent variable with ix and comparing the result with Eqs.(29a,b). The Euler's divine formula:
ei + 1 = 0 ---------- (35)
|
|
is a straight forward application of Eq.(34). It contains all the special symbols in mathematics: e, i = , , 1, and 0.
|
