Ordinary Differential Equations (1): Origins and Intuition

The World of ODEs (I): Origins and Intuition

Series Introduction: This is an 18-chapter deep dive into ordinary differential equations spanning 7 months. We'll start from scratch, understanding ODEs in the most intuitive way possible, accompanied by extensive visualizations and real-world applications from 2020 (including COVID-19 epidemic modeling). Every concept originates from practical problems and is implemented in Python, making mathematics truly come alive!

Introduction: Change Is Everywhere

If you observe the world around you, you'll notice everything is in constant flux: - Coffee is cooling - Populations are growing - Pendulums are swinging - Viruses are spreading - Hearts are beating

Differential equations are the mathematical language for describing "change." They're not just mathematical tools — they're keys to understanding the laws of nature.

Why Study Differential Equations?

In 2020, COVID-19 swept across the globe. How did scientists predict the pandemic's trajectory? How did they evaluate the effectiveness of quarantine measures? The answer lies in differential equations!

Today, we will: 1. Understand what a differential equation is (using the simplest examples) 2. Explore its historical origins 3. Examine 3 real-world applications 4. Implement our first differential equation in Python

What Is a Differential Equation? Starting with Coffee

Newton's Law of Cooling: A Story About Coffee

Imagine you've brewed a cup of coffee at 90° C and placed it in a room at 20° C. A few minutes later, you return to find the coffee has cooled.

Question: How does the coffee's temperature change over time?

Newton's Observation (1701): The rate of cooling of an object is proportional to the temperature difference between the object and its environment.

Expressed mathematically: This is a differential equation!

Interpretation: -: Temperature at time -: Rate of temperature change (derivative) -: Ambient temperature (20° C) -: Cooling constant (depends on cup material and size)

Simulating Coffee Cooling with Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Set up fonts for plotting
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False

# Define the differential equation
def coffee_cooling(T, t, k, T_env):
    """
    Newton's Law of Cooling
    T: Current temperature
    t: Time
    k: Cooling constant
    T_env: Ambient temperature
    """
    dTdt = -k * (T - T_env)
    return dTdt

# Parameter settings
T0 = 90.0        # Initial temperature (° C)
T_env = 20.0     # Ambient temperature (° C)
k = 0.1          # Cooling constant (1/minute)
t = np.linspace(0, 60, 500)  # 0-60 minutes

# Solve the differential equation
T_solution = odeint(coffee_cooling, T0, t, args=(k, T_env))

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Left plot: Temperature over time
ax1.plot(t, T_solution, 'b-', linewidth=2.5, label='Coffee Temperature')
ax1.axhline(y=T_env, color='r', linestyle='--', linewidth=2, label=f'Room Temp = {T_env}° C')
ax1.fill_between(t, T_solution.flatten(), T_env, alpha=0.3, color='orange')
ax1.set_xlabel('Time (minutes)', fontsize=12)
ax1.set_ylabel('Temperature (° C)', fontsize=12)
ax1.set_title('Newton\'s Law of Cooling', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Mark key points
ax1.plot(0, T0, 'go', markersize=10, label='Start')
ax1.text(2, T0-2, f'Initial: {T0}° C', fontsize=10, color='green')

# Right plot: Cooling rate
dT_dt = -k * (T_solution.flatten() - T_env)
ax2.plot(t, dT_dt, 'r-', linewidth=2.5)
ax2.set_xlabel('Time (minutes)', fontsize=12)
ax2.set_ylabel('Cooling Rate (° C/min)', fontsize=12)
ax2.set_title('Rate of Temperature Change', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)

plt.tight_layout()
plt.savefig('coffee_cooling.png', dpi=300, bbox_inches='tight')
plt.close()

print("✓ Figure generated: Coffee Cooling")
print(f"Initial temperature: {T0}° C")
print(f"After 10 minutes: {T_solution[np.abs(t-10).argmin()][0]:.2f}° C")
print(f"After 30 minutes: {T_solution[np.abs(t-30).argmin()][0]:.2f}° C")

Key Insights: 1. Temperature approaches ambient temperature via exponential decay 2. Cooling rate decreases over time (smaller temperature difference leads to slower cooling) 3. Theoretically, the coffee never quite reaches ambient temperature (but gets very close)

The Exact Solution

For this simple equation, we can find the exact solution:Verification: - At:✓ - As:✓

Classification and Basic Concepts

What Makes It an "Ordinary" Differential Equation?

Definition: A differential equation containing only one independent variable (usually time).

Comparison: - Ordinary Differential Equation (ODE):(one variable) - Partial Differential Equation (PDE):(multiple variables)

Order: The Highest Derivative

First-order ODE: Highest derivative is first-orderExamples: Population growth, radioactive decay, RC circuits

Second-order ODE: Highest derivative is second-orderExamples: Spring oscillation, planetary motion, RLC circuits

Linear vs. Nonlinear

Linear ODE:and its derivatives appear to the first power only

Nonlinear ODE: Contains terms like,,

Why does this matter? - Linear equations: Well-developed theory, easier to solve - Nonlinear equations: More realistic, but often require numerical methods

A Historical Journey: The Birth of Differential Equations

Newton and Leibniz (1600s)

1666: Newton invented calculus and immediately applied it to describe motion.

Newton's Second Law is itself a differential equation:This is the most important differential equation in physics!

The Bernoulli Family (1700s)

Jakob Bernoulli and Johann Bernoulli (brothers) studied many classical ODEs: - The catenary problem (shape of a hanging chain) - The brachistochrone problem (fastest descent path)

Daniel Bernoulli (their nephew/son): - Studied elastic body vibrations - Proposed the superposition principle

Euler (1700s)

1739: Euler systematized ODE theory, introducing: - General solutions for constant-coefficient linear equations - Euler's method (the first numerical algorithm!)

Laplace (1800s)

1785: The Laplace Transform - Converts differential equations to algebraic equations - An engineer's favorite tool!

Poincar é (1880s)

Revolutionary shift: Instead of seeking exact solutions, study the properties of solutions - Stability theory - Phase space - Laid the foundation for modern dynamical systems

# Visualization: Historical Timeline
fig, ax = plt.subplots(1, 1, figsize=(14, 6))

# Historical events
events = [
    (1666, "Newton\nInvents Calculus", 'blue'),
    (1690, "Bernoulli Brothers\nBrachistochrone", 'green'),
    (1739, "Euler\nSystematizes ODEs", 'red'),
    (1785, "Laplace\nTransform", 'purple'),
    (1881, "Poincar é\nQualitative Theory", 'orange'),
    (1963, "Lorenz\nChaos Theory", 'brown'),
    (2020, "COVID-19\nEpidemic Models", 'crimson')
]

years = [e[0] for e in events]
names = [e[1] for e in events]
colors = [e[2] for e in events]

# Draw timeline
ax.plot([1650, 2030], [0, 0], 'k-', linewidth=2)

for i, (year, name, color) in enumerate(events):
    # Alternate positions above and below
    y_pos = 0.3 if i % 2 == 0 else -0.3
    
    # Draw points
    ax.plot(year, 0, 'o', markersize=15, color=color, zorder=3)
    
    # Draw connecting lines
    ax.plot([year, year], [0, y_pos*0.8], color=color, linewidth=2, linestyle='--')
    
    # Add labels
    ax.text(year, y_pos, name, ha='center', va='bottom' if y_pos > 0 else 'top',
            fontsize=10, fontweight='bold', color=color,
            bbox=dict(boxstyle='round,pad=0.5', facecolor='white', edgecolor=color, linewidth=2))

ax.set_xlim(1640, 2040)
ax.set_ylim(-0.5, 0.5)
ax.set_xlabel('Year', fontsize=12, fontweight='bold')
ax.set_title('History of Differential Equations: Key Milestones', fontsize=14, fontweight='bold')
ax.axis('off')

plt.tight_layout()
plt.savefig('history_timeline.png', dpi=300, bbox_inches='tight')
plt.close()

print("✓ Figure generated: Historical Timeline")

Three Classic Application Cases

Case 1: Radioactive Decay

Background: Radioactive elements decay with a fixed "half-life."

Model:-: Number of atoms at time -: Decay constant

Solution:

Application: Carbon-14 dating (archaeology)

# Carbon-14 decay simulation
fig, ax = plt.subplots(1, 1, figsize=(10, 6))

t_half = 5730  # Carbon-14 half-life (years)
lambda_c14 = np.log(2) / t_half
t = np.linspace(0, 30000, 1000)
N = np.exp(-lambda_c14 * t)

ax.plot(t, N, 'b-', linewidth=2.5, label='C-14 Remaining')
ax.axhline(y=0.5, color='r', linestyle='--', linewidth=2, label='50% (1 half-life)')
ax.axhline(y=0.25, color='orange', linestyle='--', linewidth=2, label='25% (2 half-lives)')
ax.axvline(x=t_half, color='r', linestyle=':', alpha=0.5)
ax.axvline(x=2*t_half, color='orange', linestyle=':', alpha=0.5)

ax.set_xlabel('Time (years)', fontsize=12)
ax.set_ylabel('Fraction Remaining', fontsize=12)
ax.set_title('Radioactive Decay: Carbon-14', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

# Label half-life
ax.text(t_half, 0.5, f'  {t_half} years', fontsize=10, va='bottom', color='red')

plt.tight_layout()
plt.savefig('radioactive_decay.png', dpi=300, bbox_inches='tight')
plt.close()

print("✓ Figure generated: Radioactive Decay")

Case 2: Exponential Population Growth (Malthusian Model)

1798: Thomas Malthus proposed the simplest population model.

Assumptions: Constant birth rate, constant death rate-: Population -: Net growth rate (birth rate - death rate)

Solution:

Problem: Exponential growth is unsustainable! (Resources are limited)

Improvement: Logistic model (next chapter)

# Population growth comparison: Different growth rates
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

t = np.linspace(0, 100, 500)
P0 = 1000  # Initial population

# Left plot: Different growth rates
growth_rates = [0.01, 0.02, 0.03, 0.04]
colors = ['blue', 'green', 'orange', 'red']

for r, color in zip(growth_rates, colors):
    P = P0 * np.exp(r * t)
    ax1.plot(t, P, color=color, linewidth=2.5, label=f'r = {r} (doubling: {np.log(2)/r:.1f} yrs)')

ax1.set_xlabel('Time (years)', fontsize=12)
ax1.set_ylabel('Population', fontsize=12)
ax1.set_title('Exponential Growth: Effect of Growth Rate', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_yscale('log')

# Right plot: Semi-log coordinates
r = 0.03
P = P0 * np.exp(r * t)
ax2.semilogy(t, P, 'b-', linewidth=2.5)
ax2.set_xlabel('Time (years)', fontsize=12)
ax2.set_ylabel('Population (log scale)', fontsize=12)
ax2.set_title('Exponential Growth on Log Scale\n(Becomes a Straight Line!)', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('population_growth.png', dpi=300, bbox_inches='tight')
plt.close()

print("✓ Figure generated: Population Growth")

Case 3: Simple Harmonic Motion (Spring-Mass System)

Hooke's Law: The restoring force of a spring is proportional to displacement

Newton's Second Law:Rearranging:whereis the angular frequency.

Solution:-: Amplitude -: Initial phase

This is simple harmonic motion!

# Simple harmonic motion simulation
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

# Parameters
m = 1.0      # Mass (kg)
k = 10.0     # Spring constant (N/m)
omega = np.sqrt(k / m)
A = 1.0      # Amplitude (m)
phi = 0      # Initial phase

t = np.linspace(0, 10, 1000)
x = A * np.cos(omega * t + phi)
v = -A * omega * np.sin(omega * t + phi)

# Top plot: Displacement vs time
ax1.plot(t, x, 'b-', linewidth=2.5, label='Displacement x(t)')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axhline(y=A, color='r', linestyle='--', linewidth=1.5, alpha=0.5, label=f'Amplitude = {A}m')
ax1.axhline(y=-A, color='r', linestyle='--', linewidth=1.5, alpha=0.5)
ax1.set_xlabel('Time (seconds)', fontsize=12)
ax1.set_ylabel('Displacement (m)', fontsize=12)
ax1.set_title(f'Simple Harmonic Motion: ω = {omega:.2f} rad/s, Period = {2*np.pi/omega:.2f}s', 
             fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Bottom plot: Phase space (displacement-velocity)
ax2.plot(x, v, 'g-', linewidth=2.5)
ax2.plot(x[0], v[0], 'ro', markersize=10, label='Start')
ax2.arrow(x[50], v[50], x[51]-x[50], v[51]-v[50], 
          head_width=0.3, head_length=0.2, fc='blue', ec='blue')
ax2.set_xlabel('Displacement x (m)', fontsize=12)
ax2.set_ylabel('Velocity v (m/s)', fontsize=12)
ax2.set_title('Phase Space Portrait (Circular Orbit)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_aspect('equal')

plt.tight_layout()
plt.savefig('simple_harmonic_motion.png', dpi=300, bbox_inches='tight')
plt.close()

print("✓ Figure generated: Simple Harmonic Motion")

Geometric Interpretation: Direction Fields

What Is a Direction Field?

For a first-order ODE:At each point in theplane, the equation tells us the slope of the solution curve at that point.

Direction field: Draw these slopes (small arrows) at grid points.

Example:

def direction_field_example():
    fig, ax = plt.subplots(1, 1, figsize=(10, 8))
    
    # Define ODE
    def f(t, y):
        return y - t
    
    # Create grid
    t_vals = np.linspace(-2, 4, 20)
    y_vals = np.linspace(-2, 4, 20)
    T, Y = np.meshgrid(t_vals, y_vals)
    
    # Calculate slope at each point
    dY = f(T, Y)
    dT = np.ones_like(dY)
    
    # Normalize arrow length
    M = np.sqrt(dT**2 + dY**2)
    dT_norm = dT / M
    dY_norm = dY / M
    
    # Plot direction field
    ax.quiver(T, Y, dT_norm, dY_norm, M, cmap='viridis', alpha=0.6)
    
    # Plot several solution curves
    from scipy.integrate import odeint
    
    def ode_func(y, t):
        return y - t
    
    t_span = np.linspace(-2, 4, 500)
    initial_conditions = [-1, 0, 1, 2, 3]
    
    for y0 in initial_conditions:
        sol = odeint(ode_func, y0, t_span)
        ax.plot(t_span, sol, 'r-', linewidth=2, alpha=0.8)
    
    ax.set_xlabel('t', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title('Direction Field: dy/dt = y - t', fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-2, 4)
    ax.set_ylim(-2, 4)
    
    plt.tight_layout()
    plt.savefig('direction_field.png', dpi=300, bbox_inches='tight')
    plt.close()
    
    print("✓ Figure generated: Direction Field")

direction_field_example()

Key Insights: - The red curves are integral curves (solutions) - They flow along the arrows - Different initial values lead to different solution curves

Initial Value Problems (IVP)

Definition

Initial value problem: Given an ODE and initial conditions, find the solution

Existence and Uniqueness Theorem (Picard-Lindel ö f): Ifsatisfies certain conditions (continuous and Lipschitz continuous), then the IVP has a unique solution.

Why Are Initial Conditions Important?

Without initial conditions, a differential equation has infinitely many solutions (general solution). With initial conditions, we determine a unique solution (particular solution).

Physical interpretation: - Initial position and velocity determine an object's entire trajectory - Initial population determines future population changes

Numerical vs. Analytical Solutions

Analytical Solutions

Express the solution exactly with a formula:

Advantages: Exact, elegant
Disadvantages: Most ODEs don't have analytical solutions!

Numerical Solutions

Approximate step by step using a computer:

Advantages: Applicable to any ODE
Disadvantages: Only an approximation, requires computational resources

Reality: 99% of ODEs in engineering and science are solved numerically!

# Comparison: Analytical vs Numerical solutions
fig, ax = plt.subplots(1, 1, figsize=(10, 6))

# ODE: dy/dt = -y, y(0) = 1
# Analytical solution: y(t) = e^(-t)

t = np.linspace(0, 5, 100)
y_analytical = np.exp(-t)

# Numerical solution: Euler's method
def euler_method(f, y0, t):
    y = np.zeros(len(t))
    y[0] = y0
    for i in range(len(t) - 1):
        h = t[i+1] - t[i]
        y[i+1] = y[i] + h * f(y[i], t[i])
    return y

def f(y, t):
    return -y

y_numerical = euler_method(f, 1.0, t)

# Plotting
ax.plot(t, y_analytical, 'b-', linewidth=3, label='Analytical Solution:$y = e^{-t}$')
ax.plot(t, y_numerical, 'r--', linewidth=2, label='Numerical Solution (Euler)', alpha=0.7)
ax.plot(t, np.abs(y_analytical - y_numerical), 'g:', linewidth=2, label='Error', alpha=0.7)

ax.set_xlabel('Time t', fontsize=12)
ax.set_ylabel('y(t)', fontsize=12)
ax.set_title('Analytical vs Numerical Solution', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('numerical_vs_analytical.png', dpi=300, bbox_inches='tight')
plt.close()

print("✓ Figure generated: Numerical vs Analytical")

Looking Ahead to 2020: From Newton to COVID-19

Preview of Epidemic Modeling

In early 2020, COVID-19 emerged. The whole world was asking: - How long will the pandemic last? - Do quarantines work? - When will we reach the peak?

The answers lie in the SIR model (focus of Chapter 13):This is a system of differential equations!

Neural ODEs

2018 Breakthrough: Using neural networks to learn differential equations!

Traditional: Design, solve forNeural ODE: Let a neural network learn Applications: Time series, generative models, physics-informed modeling

Summary and Preview

Key Points from This Chapter

1. Differential equations describe change: Temperature, population, vibration...
2. ODE vs PDE: Single variable vs multiple variables
3. Order: The highest derivative's power
4. Linear vs nonlinear: Determines solution difficulty
5. Direction fields: Visualize solution flow
6. Initial value problems: Determine unique solutions
7. Numerical solutions: The main tool in practice

Why Does This Matter?

Differential equations are: - The language of physics (Newton's laws) - Tools for biology (population models) - Foundations of engineering (control systems) - Models for economics (growth theory) - Weapons for medicine (epidemic prediction)

Preview of Next Chapter

"First-Order Differential Equations"

We will learn: - Separable equations: - Linear equations: - Exact equations: - Mixing problems: Two liquids mixing - RC circuits: Capacitor charging and discharging

Each type has specific techniques!

Exercises

Basic Problems

1. Verify thatis a solution of(whereis an arbitrary constant).

2. Coffee cools from 90° C to 60° C in 10 minutes with ambient temperature 20° C. Find the cooling constant.

3. A radioactive element has a half-life of 10 years. How much of a 100-gram sample remains after 20 years?

Advanced Problems

4. Prove: All solutions ofapproach 0 (when).

5. For, sketch the direction field and predict the solution behavior.

6. Why is the exponential population model unrealistic? Propose an improvement.

Programming Problems

7. Implement Euler's method in Python to solve,, for.

8. Create an animated GIF showing solution curves with different initial values evolving in the direction field.

9. Simulate undamped pendulum motion and plot the phase space trajectory.

References

1. Boyce, W. E., & DiPrima, R. C. (2012). Elementary Differential Equations. Wiley.
2. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos. CRC Press.
3. MIT OpenCourseWare 18.03SC - Differential Equations
4. Tenenbaum, M., & Pollard, H. (1985). Ordinary Differential Equations. Dover.
5. Imperial College COVID-19 Response Team (March 2020). Impact of NPIs.

Code Resources

All code and figures for this chapter: - GitHub: [ode-series]/chapter-1 - Jupyter Notebook: Interactive version - Animated GIFs: Direction field evolution, simple harmonic motion

---

Next Chapter: First-Order Differential Equations →

---

This is Chapter 1 of "The World of Ordinary Differential Equations" series, with 18 chapters total.