Ordinary Differential Equations
Classification, Separation of Variables, and Numerical Solution
Euler · Runge–Kutta (RK2, RK4) · Higher-Order Systems · The Simple Pendulum
1. Types of Differential Equations
A differential equation (DE) is an equation that relates an unknown function to one or more of its derivatives. DEs are classified in several complementary ways; the classification determines which solution methods apply.
1.1 Order and degree
The order of a DE is the order of the highest derivative that appears. The degree is the exponent of that highest derivative (once the equation is free of radicals and fractions in the derivatives).
$$\frac{dy}{dx}=12x^{2}y\quad(\text{first order}),\qquad \frac{d^{2}y}{dx^{2}}=xy\quad(\text{second order})$$ $$\frac{d^{2}y}{dx^{2}}=2x-4\quad(\text{first degree}),\qquad \left(\frac{dy}{dx}\right)^{3}=18x-5y\quad(\text{third degree})$$1.2 Ordinary versus partial
An ordinary differential equation (ODE) involves derivatives with respect to a single independent variable, e.g. \(d^2x/dt^2=f(x)\). A partial differential equation (PDE) has two or more independent variables and contains partial derivatives:
$$\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}=xy$$These notes concern ODEs only.
1.3 Linear versus nonlinear
A linear DE is of first degree in the dependent variable and all of its derivatives: the dependent variable and its derivatives appear only to the first power and are never multiplied together or placed inside nonlinear functions. The general first-order linear form is
$$\frac{dy}{dx}+p(x)\,y=q(x),$$where \(p\) and \(q\) are functions of \(x\) (or constants). Examples:
$$\frac{dy}{dx}+y=\cos x,\qquad \frac{dy}{dx}-\frac{2y}{x}=x^{2}e^{-x},\qquad x\frac{dy}{dx}+2y=x^{2}.$$A nonlinear DE violates this: the dependent variable or a derivative is raised to a power higher than one, multiplied by another derivative, or appears inside a nonlinear function such as \(\sin y\). For example
$$\frac{d^{2}y}{dx^{2}}-\left(\frac{dy}{dx}\right)^{2}+12y=\cos x,$$and the pendulum equation of Section 4, which contains \(\sin\theta\), is nonlinear. Nonlinear equations rarely have closed-form solutions, which is a principal motivation for numerical methods.
1.4 Homogeneous versus inhomogeneous
The word homogeneous is used in two distinct senses.
(a) Homogeneous in the degree sense. A function \(f(x,y)\) is homogeneous of degree \(n\) if \(f(\lambda x,\lambda y)=\lambda^{n}f(x,y)\) for \(\lambda\neq0\). A first-order equation \(f(x,y)\,dy+g(x,y)\,dx=0\) whose \(f\) and \(g\) are homogeneous of the same degree is called homogeneous. Examples:
$$\frac{dy}{dx}=\frac{x^{3}+y^{3}}{xy^{2}+yx^{2}},\qquad \frac{dy}{dx}=\frac{3x+y}{x-y}.$$(b) Homogeneous in the forcing sense. A linear DE is homogeneous if every term contains the dependent variable (the right-hand side is zero) and inhomogeneous if there is a non-zero term on the right-hand side. Examples of inhomogeneous equations:
$$\frac{d^{2}y}{dx^{2}}+2\frac{dy}{dx}-3y=\sin x,\qquad \frac{dy}{dx}+y=x^{2},\qquad \frac{d^{2}x}{dt^{2}}+4x=e^{-t}.$$1.5 What does "solution" mean?
An analytical solution to an ODE is a function which, when substituted into the differential equation, produces a true statement. A numerical solution is instead a list of ordered pairs \((t_n,x_n)\) that gives point-wise approximations to the actual solution.
2. Solving ODEs by Separation of Variables
A first-order ODE is separable if it can be written so that all terms in the dependent variable are on one side and all terms in the independent variable are on the other:
$$\frac{dy}{dx}=g(x)\,h(y)\;\Longrightarrow\;\frac{dy}{h(y)}=g(x)\,dx\;\Longrightarrow\;\int\frac{dy}{h(y)}=\int g(x)\,dx+C.$$The constant of integration is fixed by the initial condition. We illustrate with two examples; the second requires partial fractions.
Example 1 — Radioactive decay
A radioactive isotope decays at a rate proportional to the number of nuclei present: \(dN/dt=-kN\), \(k>0\), with \(t\) in days. Initially \(N=6.0\times10^{24}\), and at \(t=10\) days \(N=6.25\times10^{22}\). (a) Find \(k\). (b) Find \(N\) at \(t=20\) days.
Solution. Separate and integrate:
$$\frac{dN}{N}=-k\,dt\;\Rightarrow\;\ln N=-kt+c\;\Rightarrow\;N=Ae^{-kt}.$$At \(t=0\), \(N=6.0\times10^{24}\), so \(A=6.0\times10^{24}\). At \(t=10\):
$$6.25\times10^{22}=6.0\times10^{24}e^{-10k}\;\Rightarrow\;\frac{1}{96}=e^{-10k}\;\Rightarrow\;k=\frac{1}{10}\ln 96\approx 0.45643\ \text{day}^{-1}.$$(b) At \(t=20\):
$$N=6.0\times10^{24}\,e^{-0.45643\times20}\approx 6.51\times10^{20}\ \text{nuclei}.$$Because \(e^{-20k}=(1/96)^{2}\), the number of nuclei falls by a factor of 96 in each 10-day interval.
Example 2 — Fox population (logistic growth)
The number of foxes \(N\) (in thousands) living in an urban area obeys \(dN/dt=2N-N^{2}\), \(t\ge0\) (years). At \(t=0\) there were 1000 foxes, i.e. \(N(0)=1\). Find \(N(t)\) and the long-term population.
Solution. Separate variables:
$$\frac{dN}{N(2-N)}=dt.$$Use partial fractions, \(\dfrac{1}{N(2-N)}=\dfrac{A}{N}+\dfrac{B}{2-N}\), so \(A(2-N)+BN=1\). Setting \(N=0\) gives \(A=\tfrac12\); setting \(N=2\) gives \(B=\tfrac12\). Hence
$$\int\left(\frac{1}{N}+\frac{1}{2-N}\right)dN=\int 2\,dt\;\Rightarrow\;\ln N-\ln(2-N)=2t+c\;\Rightarrow\;\frac{N}{2-N}=Ae^{2t}.$$At \(t=0\), \(N=1\), so \(A=1\). Solving for \(N\):
$$\frac{N}{2-N}=e^{2t}\;\Rightarrow\;N+Ne^{2t}=2e^{2t}\;\Rightarrow\;N(t)=\frac{2e^{2t}}{1+e^{2t}}=\frac{2}{1+e^{-2t}}.$$As \(t\to\infty\), \(N\to2\): the population levels off at the carrying capacity of 2000 foxes.
3. Numerical Solution of First-Order ODEs
Most ODEs cannot be solved in closed form. Numerical methods instead advance an approximate solution in small steps of size \(\Delta t\) (also written \(h\)). Throughout this section we consider the initial-value problem
$$\dot{x}=f(t,x),\qquad x(t_0)=x_0,$$and write \(t_n=t_0+n\Delta t\) and \(x_n\approx x(t_n)\). The idea in every method is the same: if \(x_n\) is known, compute slopes of the solution and use them to estimate \(x_{n+1}\).
3.1 Euler's method
Start from the finite-difference approximation of the derivative,
$$\dot{x}(t)=\frac{x(t+h)-x(t)}{h}+O(h),\qquad h=\Delta t,$$which rearranges to \(x(t+h)\approx x(t)+h\,\dot{x}(t)\). Geometrically, \(dx/dt\) is the slope of the tangent line: compute the slope at \((t_n,x_n)\) and follow it for one step.
$$x_1=x_0+\Delta t\,f(t_0,x_0),\quad x_2=x_1+\Delta t\,f(t_1,x_1),\quad\ldots,\quad x_{n+1}=x_n+\Delta t\,f(t_n,x_n).$$
Euler's method is first order: the error per step is \(O(\Delta t^{2})\) and the global error at a fixed time is \(O(\Delta t)\). It is simple but needs very small steps for good accuracy, because the slope is sampled only at the start of each interval.
3.2 Second-order Runge–Kutta (RK2)
RK2 improves on Euler by sampling the slope at more than one point in the interval and combining them. Two common versions appear in the notes.
(a) Modified Euler / Heun's method (predictor–corrector). The slope at \(t_n\) is \(f(t_n,x_n)\) and the slope at \(t_n+\Delta t\) is \(f(t_n+\Delta t,x_{n+1})\). Using the average of the two gives
$$x_{n+1}=x_n+\frac{\Delta t}{2}\Big[f(t_n,x_n)+f(t_n+\Delta t,\,x_{n+1})\Big].$$But \(x_{n+1}\) is unknown on the right-hand side, so it is first estimated with an Euler predictor step, \(x^{P}_{n+1}=x_n+\Delta t\,f(t_n,x_n)\), and then refined by the corrector step. In stage form:
$$k_1=\Delta t\,f(t_n,x_n),\quad k_2=\Delta t\,f(t_n+\Delta t,\,x_n+k_1),\quad x_{n+1}=x_n+\tfrac{1}{2}(k_1+k_2).$$(b) Midpoint method. Compute the slope at \((t_n,x_n)\), follow it for half a step to a midway point, evaluate the slope there, and use that new slope for the full step:
$$m_n=f(t_n,x_n),\quad x_{\text{mid}}=x_n+\frac{\Delta t}{2}m_n,\quad x_{n+1}=x_n+\Delta t\,f\!\left(t_n+\frac{\Delta t}{2},\,x_{\text{mid}}\right).$$Equivalently, \(k_1=f(t_n,x_n)\), \(k_2=f\!\left(t_n+\tfrac{\Delta t}{2},\,x_n+\tfrac{\Delta t}{2}k_1\right)\), \(x_{n+1}=x_n+\Delta t\,k_2\).
Both versions are second order: local error \(O(\Delta t^{3})\), global error \(O(\Delta t^{2})\). Halving the step cuts the error by about four. They cost two evaluations of \(f\) per step.
3.3 Fourth-order Runge–Kutta (RK4)
The classical RK4 method samples four slopes per step. If \(x(t_n)=x_n\) is known, then \(x(t_{n+1})\) is computed from:
| \(k_1=\Delta t\,f(t_n,x_n)\) | slope at the beginning of the interval |
| \(k_2=\Delta t\,f\!\left(t_n+\tfrac{\Delta t}{2},\,x_n+\tfrac{k_1}{2}\right)\) | slope at the midpoint using \(x_n\) and \(k_1\) |
| \(k_3=\Delta t\,f\!\left(t_n+\tfrac{\Delta t}{2},\,x_n+\tfrac{k_2}{2}\right)\) | slope at the midpoint using \(x_n\) and \(k_2\) |
| \(k_4=\Delta t\,f(t_n+\Delta t,\,x_n+k_3)\) | slope at the end of the interval using \(x_n\) and \(k_3\) |
The next value is a weighted average of the four slopes, with the two midpoint estimates counted twice:
$$x_{n+1}=x_n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right).$$RK4 is fourth order: global error \(O(\Delta t^{4})\), at a cost of four evaluations of \(f\) per step. Halving the step reduces the error by about a factor of 16.
3.4 Comparing the methods
To see the orders of accuracy in practice, we solve \(\dot{x}=x\), \(x(0)=1\) (exact solution \(e^{t}\)) and report the absolute error at \(t=1\) for several step sizes.
| \(\Delta t\) | Euler (order 1) | RK2 (order 2) | RK4 (order 4) |
|---|---|---|---|
| 0.5 | 4.68e-01 | 7.77e-02 | 9.36e-04 |
| 0.25 | 2.77e-01 | 2.34e-02 | 7.19e-05 |
| 0.1 | 1.25e-01 | 4.20e-03 | 2.08e-06 |
| 0.05 | 6.50e-02 | 1.09e-03 | 1.36e-07 |
Table 1. Absolute error at \(t=1\) for \(\dot{x}=x\), \(x(0)=1\).
The errors shrink roughly in proportion to \(\Delta t\), \(\Delta t^{2}\) and \(\Delta t^{4}\) respectively. RK4 is usually the best trade-off between accuracy and effort, which is why it is the standard general-purpose fixed-step method.
4. Higher-Order ODEs and the Simple Pendulum
The methods of Section 3 apply to first-order equations. A higher-order ODE is handled by rewriting it as a system of first-order ODEs, which the same methods then solve component by component.
4.1 Reduction to a first-order system
For a second-order ODE \(\ddot{x}=f(t,x,\dot{x})\) introduce a new variable \(u=\dot{x}\). Then
$$\dot{x}=u,\qquad \dot{u}=f(t,x,u).$$For a third-order ODE \(d^{3}x/dt^{3}=f(t,x,\dot{x},\ddot{x})\) set \(u=\dot{x}\) and \(v=\ddot{x}\):
$$\dot{x}=u,\qquad \dot{u}=v,\qquad \dot{v}=f(t,x,u,v).$$In general an order-\(m\) ODE becomes \(m\) coupled first-order equations. The state is a vector, and the numerical method is applied to the vector.
4.2 Methods for systems
For a pair \(\dot{x}=f(t,x,y)\), \(\dot{y}=g(t,x,y)\), the slopes for both variables must be computed at the same stage before either variable is advanced. With the modified Euler (RK2) method:
$$\begin{aligned}k_1&=\Delta t\,f(t_n,x_n,y_n), & l_1&=\Delta t\,g(t_n,x_n,y_n),\\ k_2&=\Delta t\,f(t_n+\Delta t,x_n+k_1,y_n+l_1), & l_2&=\Delta t\,g(t_n+\Delta t,x_n+k_1,y_n+l_1),\end{aligned}$$ $$x_{n+1}=x_n+\tfrac12(k_1+k_2),\qquad y_{n+1}=y_n+\tfrac12(l_1+l_2).$$RK4 extends in exactly the same way, with four stages for each variable.
4.3 The simple pendulum
Consider a pendulum of length \(L\) with a bob swinging through angle \(\theta\) (in radians) from the vertical. The arc length is \(x=L\theta\). Only the component of gravity along the arc, \(-g\sin\theta\), accelerates the bob, so
$$a=\ddot{x}=L\ddot{\theta}=-g\sin\theta\;\;\Longrightarrow\;\;\ddot{\theta}=-\frac{g}{L}\sin\theta.$$For small angles, \(\sin\theta\approx\theta\) and the motion is harmonic with period \(T=2\pi\sqrt{L/g}\). For large angles the equation is nonlinear and has no elementary closed-form solution, so we solve it numerically. Air resistance, proportional to the angular velocity with coefficient \(\mu\), adds a damping term:
$$\ddot{\theta}=-\mu\dot{\theta}-\frac{g}{L}\sin\theta.$$Phase space. The state of the pendulum at any instant is the vector \((\theta,\dot\theta)\). Plotting \(\dot\theta\) against \(\theta\) gives the phase space picture, in which the solution traces a curve. Writing \(\omega=\dot\theta\) breaks the single second-order ODE into two first-order ODEs:
$$\frac{d}{dt}\begin{bmatrix}\theta\\ \omega\end{bmatrix}=\begin{bmatrix}\omega\\ -\mu\omega-\dfrac{g}{L}\sin\theta\end{bmatrix}$$ $$\frac{d\theta}{dt}=\omega=f_1(t,\theta,\omega),\qquad \frac{d\omega}{dt}=-\mu\omega-\frac{g}{L}\sin\theta=f_2(t,\theta,\omega).$$Specify initial conditions, for example \(\theta_0=\pi/3\) and \(\omega_0=0\) (released from rest).
RK4 applied to the pendulum. Let \(h=\Delta t\). If \(\theta_n\) and \(\omega_n\) are known at \(t_n\), compute
$$\begin{aligned} k_{1\theta}&=f_1(t_n,\theta_n,\omega_n), & k_{1\omega}&=f_2(t_n,\theta_n,\omega_n),\\ k_{2\theta}&=f_1\!\left(t_n+\tfrac h2,\ \theta_n+\tfrac h2k_{1\theta},\ \omega_n+\tfrac h2k_{1\omega}\right), & k_{2\omega}&=f_2\!\left(t_n+\tfrac h2,\ \theta_n+\tfrac h2k_{1\theta},\ \omega_n+\tfrac h2k_{1\omega}\right),\\ k_{3\theta}&=f_1\!\left(t_n+\tfrac h2,\ \theta_n+\tfrac h2k_{2\theta},\ \omega_n+\tfrac h2k_{2\omega}\right), & k_{3\omega}&=f_2\!\left(t_n+\tfrac h2,\ \theta_n+\tfrac h2k_{2\theta},\ \omega_n+\tfrac h2k_{2\omega}\right),\\ k_{4\theta}&=f_1\!\left(t_n+h,\ \theta_n+hk_{3\theta},\ \omega_n+hk_{3\omega}\right), & k_{4\omega}&=f_2\!\left(t_n+h,\ \theta_n+hk_{3\theta},\ \omega_n+hk_{3\omega}\right), \end{aligned}$$and then advance both variables with the RK4 weights:
$$\theta_{n+1}=\theta_n+\frac h6\left[k_{1\theta}+2k_{2\theta}+2k_{3\theta}+k_{4\theta}\right],\qquad \omega_{n+1}=\omega_n+\frac h6\left[k_{1\omega}+2k_{2\omega}+2k_{3\omega}+k_{4\omega}\right].$$4.4 Results
Figure 2 shows the numerical solution for \(L=1\) m, \(g=9.81\) m/s², \(\theta_0=\pi/3\), \(\omega_0=0\) and step size \(h=0.05\) s for RK4 (0.01 s for Euler). RK4 is indistinguishable from a fine-step reference solution even with the larger step, while Euler's method visibly drifts: the amplitude grows artificially, so the pendulum gains energy it should not have. In the phase plane the undamped pendulum traces a closed orbit, whereas with air resistance (\(\mu=0.5\)) the orbit spirals into the stable equilibrium \(\theta=\omega=0\).
5. Summary
- DEs are classified by order, degree, ordinary vs. partial, linear vs. nonlinear and homogeneous vs. inhomogeneous.
- Separable first-order ODEs are solved by moving all \(y\)-terms to one side, all \(x\)- (or \(t\)-) terms to the other, and integrating; initial conditions fix the constant.
- Euler follows one slope per step (first order); RK2 averages two slopes (second order); RK4 averages four slopes (fourth order).
- A higher-order ODE is converted to a system of first-order ODEs, which Euler, RK2 or RK4 can then advance as a vector. The pendulum equation \(\ddot\theta=-\mu\dot\theta-\tfrac gL\sin\theta\) becomes \(d\theta/dt=\omega\), \(d\omega/dt=-\mu\omega-\tfrac gL\sin\theta\).