Ordinary Differential Equations (2): First-Order Equations

First-order differential equations are the most fundamental and widely used type of differential equations. From compound interest calculations to drug metabolism, from chemical reactions to circuit analysis, first-order ODEs are everywhere. Mastering the techniques for solving these equations forms the foundation for understanding more complex systems.

Separable Equations: The Most Natural Form

What Is a Separable Equation?

Consider the general form of a first-order ODE:

Ifcan be written as the product of a function ofand a function of:This is a separable equation.

Core idea: Move all terms containingto the left side and all terms containingto the right side, then integrate both sides separately.

Classic Example: Exponential Growth and Decay

Example 1: Solve (whereis a constant)

Solution: This is a separable equation withand.

Separating variables:Integrating both sides: Let (including the case), we get the general solution:

Physical interpretation: -: Exponential growth (population growth, bacterial reproduction, compound interest) -: Exponential decay (radioactive decay, drug metabolism, coffee cooling)

import numpy as np
import matplotlib.pyplot as plt

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

# Exponential functions with different k values
x = np.linspace(0, 5, 200)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Exponential growth
for k in [0.5, 1.0, 1.5, 2.0]:
    axes[0].plot(x, np.exp(k * x), label=f'k = {k}', linewidth=2)
axes[0].set_xlabel('x', fontsize=12)
axes[0].set_ylabel('y', fontsize=12)
axes[0].set_title('Exponential Growth (k > 0)', fontsize=14, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0, 30)

# Exponential decay
for k in [-0.5, -1.0, -1.5, -2.0]:
    axes[1].plot(x, np.exp(k * x), label=f'k = {k}', linewidth=2)
axes[1].set_xlabel('x', fontsize=12)
axes[1].set_ylabel('y', fontsize=12)
axes[1].set_title('Exponential Decay (k < 0)', fontsize=14, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fig1_exponential.png', dpi=150, bbox_inches='tight')
plt.close()

Compound Interest Problem

Scenario: You depositdollars in a bank with annual interest rate (continuous compounding). Find the principal afteryears.

Modeling: Letbe the principal at time. The rate of change equals the interest:Initial condition: Solution:

Example calculation: Deposit$10,000 with 5% annual interest rate

def compound_interest():
    P0 = 10000  # Initial principal
    r = 0.05    # 5% annual interest rate
    t = np.linspace(0, 30, 300)
    
    # Continuous compounding
    P_continuous = P0 * np.exp(r * t)
    
    # Yearly compounding comparison
    P_yearly = P0 * (1 + r) ** t
    
    # Monthly compounding
    P_monthly = P0 * (1 + r/12) ** (12 * t)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(t, P_continuous, 'b-', linewidth=2.5, label='Continuous Compounding')
    ax.plot(t, P_monthly, 'g--', linewidth=2, label='Monthly Compounding')
    ax.plot(t, P_yearly, 'r:', linewidth=2, label='Yearly Compounding')
    
    ax.set_xlabel('Years', fontsize=12)
    ax.set_ylabel('Principal ($)', fontsize=12)
    ax.set_title(f'Compound Interest Comparison (r = {r*100}%)', fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    
    # Annotate values at 10 and 30 years
    for t_mark in [10, 30]:
        idx = np.abs(t - t_mark).argmin()
        ax.annotate(f'${P_continuous[idx]:.0f}', 
                   xy=(t_mark, P_continuous[idx]),
                   xytext=(t_mark+2, P_continuous[idx]+2000),
                   fontsize=10, arrowprops=dict(arrowstyle='->', color='blue'))
    
    plt.tight_layout()
    plt.savefig('fig2_compound_interest.png', dpi=150, bbox_inches='tight')
    plt.close()

compound_interest()

Key insight:

Time	Yearly Compounding	Continuous Compounding	Difference
10 years	$16,288.95	$16,487.21	$198.26
30 years	$43,219.42	$44,816.89	$1,597.47

Continuous compounding seems to yield only slightly more interest, but over the long term, the difference becomes significant. This is the magic of

Drug Metabolism: The Half-Life Concept

Background: You take a pill, and your body metabolizes the drug at a constant rate. How long until the effect is halved?

Model: Letbe the drug concentration in blood

Solution: Half-life: The time for concentration to drop to half

Example: Ibuprofen has a half-life of about 2 hours. If you initially take 400mg:

def drug_metabolism():
    C0 = 400  # Initial concentration (mg)
    t_half = 2  # Half-life (hours)
    k = np.log(2) / t_half
    
    t = np.linspace(0, 12, 200)
    C = C0 * np.exp(-k * t)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(t, C, 'b-', linewidth=2.5, label='Drug Concentration')
    ax.axhline(y=200, color='r', linestyle='--', label='50% (1 half-life)')
    ax.axhline(y=100, color='orange', linestyle='--', label='25% (2 half-lives)')
    ax.axhline(y=50, color='green', linestyle='--', label='12.5% (3 half-lives)')
    
    # Mark half-lives
    for i in range(1, 4):
        t_i = i * t_half
        C_i = C0 / (2**i)
        ax.plot(t_i, C_i, 'ko', markersize=8)
        ax.text(t_i + 0.2, C_i + 15, f't = {t_i}h\nC = {C_i}mg', fontsize=9)
    
    ax.set_xlabel('Time (hours)', fontsize=12)
    ax.set_ylabel('Drug Concentration (mg)', fontsize=12)
    ax.set_title('Drug Metabolism: Ibuprofen (Half-life = 2h)', fontsize=14, fontweight='bold')
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 450)
    
    plt.tight_layout()
    plt.savefig('fig3_drug_metabolism.png', dpi=150, bbox_inches='tight')
    plt.close()

drug_metabolism()

Clinical significance: - Typically after 3-5 half-lives, the drug is almost completely metabolized - Dosing intervals are usually designed based on half-life - Ibuprofen: Take every 4-6 hours (about 2-3 half-lives)

Logistic Equation: Growth Under Resource Constraints

Problem with Malthusian model: Exponential growthleads to unlimited population growth, which is clearly unrealistic.

Improvement: Pierre Fran ç ois Verhulst (1838) proposed the Logistic model -: Intrinsic growth rate -: Carrying capacity

Intuition: - When:(exponential growth) - When:(growth stops) - When:(population declines)

Solution: This is a separable equationUsing partial fraction decomposition:Integrating yields: Rearranging gives the Logistic function:

def logistic_growth():
    K = 1000  # Carrying capacity
    r = 0.3   # Growth rate
    
    t = np.linspace(0, 30, 300)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Logistic curves with different initial values
    P0_values = [10, 50, 100, 500, 1200]
    colors = plt.cm.viridis(np.linspace(0, 1, len(P0_values)))
    
    for P0, color in zip(P0_values, colors):
        P = K / (1 + (K/P0 - 1) * np.exp(-r * t))
        axes[0].plot(t, P, color=color, linewidth=2, label=f'P ₀ = {P0}')
    
    axes[0].axhline(y=K, color='red', linestyle='--', linewidth=2, label=f'K = {K}')
    axes[0].set_xlabel('Time', fontsize=12)
    axes[0].set_ylabel('Population', fontsize=12)
    axes[0].set_title('Logistic Growth: Different Initial Conditions', fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # Growth rate vs population
    P_range = np.linspace(0, 1.5 * K, 200)
    dP_dt = r * P_range * (1 - P_range / K)
    
    axes[1].plot(P_range, dP_dt, 'b-', linewidth=2.5)
    axes[1].axvline(x=K/2, color='green', linestyle='--', label='Maximum growth at P = K/2')
    axes[1].axvline(x=K, color='red', linestyle='--', label='Carrying capacity K')
    axes[1].axhline(y=0, color='black', linewidth=0.5)
    axes[1].fill_between(P_range[P_range <= K], dP_dt[P_range <= K], alpha=0.3, color='blue')
    axes[1].fill_between(P_range[P_range > K], dP_dt[P_range > K], alpha=0.3, color='red')
    
    axes[1].set_xlabel('Population P', fontsize=12)
    axes[1].set_ylabel('Growth Rate dP/dt', fontsize=12)
    axes[1].set_title('Growth Rate vs Population', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig4_logistic.png', dpi=150, bbox_inches='tight')
    plt.close()

logistic_growth()

Characteristics of the S-shaped curve:

1. Inflection point: Maximum growth rate at

2. Asymptote:(never exceeds by much)

3. Symmetry: Symmetric about the inflection point

Real-world applications: - Bacterial growth in a petri dish - Technology adoption curves - Epidemic spread (early stages) - Market saturation models

First-Order Linear Equations: Integrating Factor Method

Standard Form

The standard form of a first-order linear ODE:

Characteristic: Bothandappear to the first power, and coefficients can be functions of.

Motivation for the Integrating Factor Method

Consider the derivative of a product:If we can find a functionsuch that:Then the equation becomes easy to integrate!

Condition: Comparing coefficients, we needThis is itself a separable equation: This is the integrating factor!

General Solution Formula

Multiply both sides of the original equation by:Integrate:

General solution:

Detailed Example

Example 2: Solve Identify:, Calculate integrating factor:

Multiply by integrating factor:The left side is:

Integrate: (using substitution)

General solution:

def linear_ode_example():
    x = np.linspace(-3, 3, 300)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # Solution curves for different C values
    C_values = [-2, -1, -0.5, 0, 0.5, 1, 2]
    colors = plt.cm.coolwarm(np.linspace(0, 1, len(C_values)))
    
    for C, color in zip(C_values, colors):
        y = 0.5 + C * np.exp(-x**2)
        ax.plot(x, y, color=color, linewidth=2, label=f'C = {C}')
    
    # Equilibrium solution
    ax.axhline(y=0.5, color='black', linestyle='--', linewidth=2, label='Equilibrium y = 1/2')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title("Solution Family: y' + 2xy = x", fontsize=14, fontweight='bold')
    ax.legend(loc='upper right')
    ax.grid(True, alpha=0.3)
    ax.set_ylim(-2, 3)
    
    plt.tight_layout()
    plt.savefig('fig5_linear_ode.png', dpi=150, bbox_inches='tight')
    plt.close()

linear_ode_example()

Observation: All solution curves approach the equilibrium solutionThis is becauseas.

RC Circuit: A Classic Application of First-Order Linear ODEs

Scenario: A resistorand capacitorin series, connected to a voltage source. Find the voltage across the capacitor.

Circuit equation (Kirchhoff's voltage law):Since:Rearranging to standard form:Define the time constant:

Case 1: Constant voltage source

Integrating factor:General solution:

Physical interpretation: - Capacitor voltage exponentially approachesfrom initial voltage - Time constantdetermines charging/discharging speed - At, about 63.2% of charging is complete - At, almost fully charged (99.3%)

def rc_circuit():
    R = 1000  # 1k Ω
    C = 0.001  # 1mF
    tau = R * C  # 1s
    V0 = 5  # Source voltage 5V
    
    t = np.linspace(0, 5 * tau, 300)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Charging process (initial 0V)
    Vc_charge = V0 * (1 - np.exp(-t / tau))
    i_charge = (V0 / R) * np.exp(-t / tau)
    
    axes[0].plot(t, Vc_charge, 'b-', linewidth=2.5, label='Capacitor Voltage$V_C$')
    axes[0].plot(t, V0 - Vc_charge, 'r--', linewidth=2, label='Resistor Voltage$V_R$')
    axes[0].axhline(y=V0, color='green', linestyle=':', label=f'Source Voltage = {V0}V')
    axes[0].axvline(x=tau, color='gray', linestyle='--', alpha=0.5)
    axes[0].text(tau, 0.2, f'τ = {tau}s', fontsize=10)
    
    # Mark 63.2% point
    axes[0].plot(tau, V0 * 0.632, 'ko', markersize=8)
    axes[0].annotate('63.2%', xy=(tau, V0*0.632), xytext=(tau+0.5, V0*0.5),
                    fontsize=10, arrowprops=dict(arrowstyle='->', color='black'))
    
    axes[0].set_xlabel('Time (s)', fontsize=12)
    axes[0].set_ylabel('Voltage (V)', fontsize=12)
    axes[0].set_title('RC Circuit: Charging', fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # Discharging process (initial 5V, disconnected from source)
    Vc_discharge = V0 * np.exp(-t / tau)
    
    axes[1].plot(t, Vc_discharge, 'b-', linewidth=2.5, label='Capacitor Voltage$V_C$')
    axes[1].axhline(y=0, color='green', linestyle=':', label='Ground')
    axes[1].axvline(x=tau, color='gray', linestyle='--', alpha=0.5)
    
    # Mark 36.8% point
    axes[1].plot(tau, V0 * 0.368, 'ko', markersize=8)
    axes[1].annotate('36.8%', xy=(tau, V0*0.368), xytext=(tau+0.5, V0*0.5),
                    fontsize=10, arrowprops=dict(arrowstyle='->', color='black'))
    
    axes[1].set_xlabel('Time (s)', fontsize=12)
    axes[1].set_ylabel('Voltage (V)', fontsize=12)
    axes[1].set_title('RC Circuit: Discharging', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig6_rc_circuit.png', dpi=150, bbox_inches='tight')
    plt.close()

rc_circuit()

Exact Equations and Integrating Factors

What Is an Exact Equation?

Consider an equation of the formOr equivalently in differential form:If there exists a functionsuch that:That is,and, then the equation is exact.

Test for exactness:This is because(equality of mixed partial derivatives).

Solving Exact Equations

If the equation is exact, the general solution is.

Steps: 1. Verify$ = = MxF = M , dx + g(y) = Ng(y)F(x,y) = C$ Example 3: Solve Verify exactness: -, -,They're equal! The equation is exact.

Find:Using:So, thus(constant can be ignored).

General solution:

Making Equations Exact: Integrating Factors

When an equation is not exact (), we can try to find an integrating factorsuch that:becomes an exact equation.

Special case 1:(depends only on)

Requiresto depend only on.

Then:

Special case 2:(depends only on)

Requiresto depend only on.

Then:

Bernoulli Equations

Form and Reduction

Bernoulli equation:where(otherwise it's a linear equation).

Technique: Make the substitutionThenDivide the original equation by:That is:This is a linear equation in Example 4: SolveThis is a Bernoulli equation with.

Let, thenDivide the original equation by: This is a linear equation! Integrating factor Using integration by parts: Substituting back:

Mixing Problems: Salt in a Tank

Problem Setup

A 1000-liter tank initially contains pure water. Salt water with concentration 2 g/L flows in at 5 L/min, and the mixed solution flows out at 5 L/min. Find the amount of salt at any time.

Modeling

Letbe the grams of salt in the tank at time.

Inflow rate of salt:g/min

Outflow rate of salt:g/min

(because outflow concentration = current tank concentration =g/L)

Differential equation:

Initial condition:(pure water)

Solution

This is a first-order linear equation. Rearranging to standard form:Integrating factor: Using:, so Solution:

def salt_tank():
    t = np.linspace(0, 1500, 300)
    Q = 2000 * (1 - np.exp(-t / 200))
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Salt amount over time
    axes[0].plot(t, Q, 'b-', linewidth=2.5, label='Salt Amount Q(t)')
    axes[0].axhline(y=2000, color='red', linestyle='--', linewidth=2, label='Equilibrium = 2000g')
    
    # Mark key moments
    for t_mark in [200, 400, 600, 1000]:
        Q_mark = 2000 * (1 - np.exp(-t_mark / 200))
        axes[0].plot(t_mark, Q_mark, 'ko', markersize=6)
        axes[0].annotate(f't={t_mark}min\nQ={Q_mark:.0f}g', 
                        xy=(t_mark, Q_mark), xytext=(t_mark+50, Q_mark-200),
                        fontsize=9)
    
    axes[0].set_xlabel('Time (minutes)', fontsize=12)
    axes[0].set_ylabel('Salt Amount (grams)', fontsize=12)
    axes[0].set_title('Salt Tank Problem', fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # Concentration change
    C = Q / 1000  # g/L
    axes[1].plot(t, C, 'g-', linewidth=2.5, label='Concentration C(t)')
    axes[1].axhline(y=2, color='red', linestyle='--', linewidth=2, label='Inlet Concentration = 2 g/L')
    
    axes[1].set_xlabel('Time (minutes)', fontsize=12)
    axes[1].set_ylabel('Concentration (g/L)', fontsize=12)
    axes[1].set_title('Salt Concentration Over Time', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig7_salt_tank.png', dpi=150, bbox_inches='tight')
    plt.close()

salt_tank()

Key insights: - Salt amount approaches 2000g (concentration approaches 2 g/L, equal to inlet concentration) - Time constantminutes - After about 1000 minutes (about 17 hours), equilibrium is essentially reached

Variation: Changing Tank Volume

If inflow rateoutflow rate, the tank volume changes!

Example: Inflow 6 L/min (2 g/L), outflow 4 L/min, initial 500L pure water.

Volume:liters

Differential equation:This requires a more complex integrating factor, but the principle is the same.

Existence and Uniqueness Theorem

Picard-Lindel ö f Theorem

Theorem: Consider the initial value problemIfsatisfies the following in a rectangular regioncontaining: 1.is continuous with respect to$x, yy$(Lipschitz condition)

Then a unique solution exists near.

Why Is This Important?

Positive side: Tells us when we can safely solve the equation

Negative side: When conditions are not met, we may have: - No solution - Infinitely many solutions - Solution exists only on a finite interval

Example 5:,Check: is unbounded atThis initial value problem has infinitely many solutions: -(trivial solution) -for any(non-trivial solutions) - - Piecewise functions...

def non_unique_solution():
    x = np.linspace(-1, 3, 300)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # y = 0 (trivial solution)
    ax.axhline(y=0, color='blue', linewidth=2.5, label='y = 0 (trivial)')
    
    # y = x^3
    ax.plot(x, x**3, 'r-', linewidth=2, label='y = x ³')
    
    # y = (x-1)^3 for x >= 1, 0 otherwise
    y1 = np.where(x >= 1, (x-1)**3, 0)
    ax.plot(x, y1, 'g--', linewidth=2, label='y = (x-1)³ for x ≥ 1')
    
    # y = (x-2)^3 for x >= 2, 0 otherwise
    y2 = np.where(x >= 2, (x-2)**3, 0)
    ax.plot(x, y2, 'm:', linewidth=2, label='y = (x-2)³ for x ≥ 2')
    
    ax.plot(0, 0, 'ko', markersize=10, label='Initial point (0, 0)')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title("Non-uniqueness: dy/dx = 3y^(2/3), y(0) = 0", fontsize=14, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-1, 3)
    ax.set_ylim(-2, 10)
    
    plt.tight_layout()
    plt.savefig('fig8_non_unique.png', dpi=150, bbox_inches='tight')
    plt.close()

non_unique_solution()

Introduction to Numerical Methods

Euler's Method

Idea: Approximate the curve with tangent lineswhereis the step size.

Geometric interpretation: Starting from, take a small step in the direction of slope.

def euler_method_demo():
    # ODE: dy/dx = y, y(0) = 1
    # Exact solution: y = e^x
    
    def f(x, y):
        return y
    
    def euler(f, x0, y0, h, n_steps):
        x = [x0]
        y = [y0]
        for _ in range(n_steps):
            y_new = y[-1] + h * f(x[-1], y[-1])
            x_new = x[-1] + h
            x.append(x_new)
            y.append(y_new)
        return np.array(x), np.array(y)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    x_exact = np.linspace(0, 2, 200)
    y_exact = np.exp(x_exact)
    
    # Compare different step sizes
    step_sizes = [0.5, 0.2, 0.1, 0.05]
    colors = ['red', 'orange', 'green', 'purple']
    
    for h, color in zip(step_sizes, colors):
        n_steps = int(2 / h)
        x_euler, y_euler = euler(f, 0, 1, h, n_steps)
        axes[0].plot(x_euler, y_euler, 'o-', color=color, markersize=4, 
                    linewidth=1.5, label=f'h = {h}')
    
    axes[0].plot(x_exact, y_exact, 'b-', linewidth=2.5, label='Exact:$y = e^x$')
    axes[0].set_xlabel('x', fontsize=12)
    axes[0].set_ylabel('y', fontsize=12)
    axes[0].set_title("Euler's Method: Different Step Sizes", fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # Error analysis
    step_sizes_log = np.logspace(-2, 0, 20)
    errors = []
    
    for h in step_sizes_log:
        n_steps = int(2 / h)
        x_euler, y_euler = euler(f, 0, 1, h, n_steps)
        error = np.abs(y_euler[-1] - np.exp(2))
        errors.append(error)
    
    axes[1].loglog(step_sizes_log, errors, 'bo-', linewidth=2, markersize=6, label='Error')
    axes[1].loglog(step_sizes_log, step_sizes_log, 'r--', linewidth=2, label='O(h) reference')
    axes[1].set_xlabel('Step Size h', fontsize=12)
    axes[1].set_ylabel('Error at x = 2', fontsize=12)
    axes[1].set_title('Error vs Step Size (Log-Log)', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig('fig9_euler.png', dpi=150, bbox_inches='tight')
    plt.close()

euler_method_demo()

Error analysis: - Local truncation error:(per step) - Global error:(cumulative)

Euler's method is a first-order method.

Improved Euler Method (Heun's Method)

Idea: Use the average of slopes$

$ This is a second-order method with global error.

Fourth-Order Runge-Kutta Method (RK4)

The most commonly used numerical method:$

$

Global error:

def rk4_comparison():
    def f(x, y):
        return y
    
    def euler(f, x0, y0, h, n_steps):
        x, y = [x0], [y0]
        for _ in range(n_steps):
            y.append(y[-1] + h * f(x[-1], y[-1]))
            x.append(x[-1] + h)
        return np.array(x), np.array(y)
    
    def rk4(f, x0, y0, h, n_steps):
        x, y = [x0], [y0]
        for _ in range(n_steps):
            k1 = f(x[-1], y[-1])
            k2 = f(x[-1] + h/2, y[-1] + h*k1/2)
            k3 = f(x[-1] + h/2, y[-1] + h*k2/2)
            k4 = f(x[-1] + h, y[-1] + h*k3)
            y.append(y[-1] + h/6 * (k1 + 2*k2 + 2*k3 + k4))
            x.append(x[-1] + h)
        return np.array(x), np.array(y)
    
    h = 0.2
    n_steps = 10
    
    x_euler, y_euler = euler(f, 0, 1, h, n_steps)
    x_rk4, y_rk4 = rk4(f, 0, 1, h, n_steps)
    x_exact = np.linspace(0, 2, 200)
    y_exact = np.exp(x_exact)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    ax.plot(x_exact, y_exact, 'b-', linewidth=2.5, label='Exact')
    ax.plot(x_euler, y_euler, 'r--o', linewidth=2, markersize=6, label='Euler')
    ax.plot(x_rk4, y_rk4, 'g-s', linewidth=2, markersize=6, label='RK4')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title(f'Euler vs RK4 (h = {h})', fontsize=14, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # Display errors
    euler_error = np.abs(y_euler[-1] - np.exp(2))
    rk4_error = np.abs(y_rk4[-1] - np.exp(2))
    ax.text(0.5, 5, f'Euler error: {euler_error:.4f}\nRK4 error: {rk4_error:.6f}', 
           fontsize=11, bbox=dict(boxstyle='round', facecolor='wheat'))
    
    plt.tight_layout()
    plt.savefig('fig10_rk4.png', dpi=150, bbox_inches='tight')
    plt.close()

rk4_comparison()

Python Practice: Comprehensive Applications

Solving ODEs with SciPy

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

# Define ODE: Cooling law
def cooling(t, T, k, T_env):
    return -k * (T - T_env)

# Parameters
k = 0.1
T_env = 20
T0 = 90

# Solve
sol = solve_ivp(
    fun=lambda t, T: cooling(t, T, k, T_env),
    t_span=[0, 60],
    y0=[T0],
    method='RK45',
    dense_output=True
)

# Plot
t = np.linspace(0, 60, 300)
T = sol.sol(t)[0]

plt.figure(figsize=(10, 6))
plt.plot(t, T, 'b-', linewidth=2.5, label='Numerical Solution')
plt.plot(t, T_env + (T0 - T_env) * np.exp(-k * t), 'r--', 
         linewidth=2, label='Analytical Solution')
plt.xlabel('Time (min)', fontsize=12)
plt.ylabel('Temperature (° C)', fontsize=12)
plt.title('Cooling Law: SciPy solve_ivp', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('fig11_scipy.png', dpi=150, bbox_inches='tight')
plt.close()

Direction Field Visualization

def plot_direction_field(f, x_range, y_range, ax=None, density=20):
    """Plot direction field for a first-order ODE"""
    if ax is None:
        fig, ax = plt.subplots(figsize=(10, 8))
    
    x = np.linspace(*x_range, density)
    y = np.linspace(*y_range, density)
    X, Y = np.meshgrid(x, y)
    
    # Calculate slopes
    dY = f(X, Y)
    dX = np.ones_like(dY)
    
    # Normalize
    N = np.sqrt(dX**2 + dY**2)
    dX, dY = dX/N, dY/N
    
    ax.quiver(X, Y, dX, dY, N, cmap='coolwarm', alpha=0.6)
    return ax

# Example: dy/dx = sin(x) - y
def f(x, y):
    return np.sin(x) - y

fig, ax = plt.subplots(figsize=(10, 8))
plot_direction_field(f, (-2*np.pi, 2*np.pi), (-3, 3), ax)

# Plot several solution curves
from scipy.integrate import odeint

x_span = np.linspace(-2*np.pi, 2*np.pi, 300)
for y0 in [-2, -1, 0, 1, 2]:
    sol = odeint(lambda y, x: np.sin(x) - y, y0, x_span)
    ax.plot(x_span, sol, 'k-', linewidth=1.5, alpha=0.7)

ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title("Direction Field: dy/dx = sin(x) - y", fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('fig12_direction_field.png', dpi=150, bbox_inches='tight')
plt.close()

Summary

In this chapter, we systematically studied the main types of first-order differential equations and their solution methods:

Type	Standard Form	Solution Method
Separable		Separate variables and integrate
First-order linear		Integrating factor
Exact	,	Find potential function
Bernoulli		Substitution

Application areas: - Compound interest and investments (exponential growth) - Drug metabolism (exponential decay) - Population dynamics (Logistic equation) - Circuit analysis (RC circuits) - Chemical reactions (mixing problems)

Numerical methods: - Euler's method: Simple but low accuracy - RK4 method: Most commonly used in engineering - SciPy: Modern Python tools

Exercises

Basic Problems

1. Solve the separable equation,.

2. Solve the linear equation.

3. Determine ifis exact. If not, find an integrating factor.

4. A pool initially has 1000 liters of water containing 50kg of salt. Pure water flows in at 10 L/min, mixed solution flows out at 10 L/min. Find.

5. In an RC circuit with,, initial capacitor voltage is 0, connected to a 10V source. Find.

Advanced Problems

6. Solve the Bernoulli equation.

7. Prove: For equations of the form, lettingconverts it to a separable equation.

8. Discuss the existence and uniqueness of solutions for,.

9. Bacteria count satisfies the Logistic equation,. Find:

   • Expression for - When is?
   • Time of maximum growth

10. Solve the equation(Hint: use polar coordinates)

Programming Problems

11. Implement the improved Euler method (Heun's method) in Python and compare accuracy with Euler's method.

12. Plot the direction field for the Logistic equation with solution curves for different initial values.

13. Simulate continuous medication: Take 500mg every 6 hours, half-life is 4 hours. Plot the blood drug concentration curve.

14. Implement adaptive step-size RK45 and compare efficiency with fixed step-size methods.

References

1. Boyce, W. E., & DiPrima, R. C. (2012). Elementary Differential Equations and Boundary Value Problems. Wiley.
2. Zill, D. G. (2017). A First Course in Differential Equations with Modeling Applications. Cengage.
3. Edwards, C. H., & Penney, D. E. (2014). Differential Equations and Boundary Value Problems. Pearson.
4. SciPy Documentation: scipy.integrate

---

Next Chapter Preview: Higher-Order Linear Differential Equations — Spring oscillations, RLC circuits, resonance phenomena