Introduction to Differential Equations

Nearly every law of nature is written as a differential equation. Newton's second law \(F = ma\) becomes \(m\ddot{x} = F(x, \dot{x}, t)\) the moment you express acceleration as the second derivative of position. The heat equation, Maxwell's equations, the Black-Scholes model, population dynamics, neural network training — all are differential equations in disguise.

Before you can solve any of these, you need a precise vocabulary: What is a differential equation? What does "order" mean? What distinguishes a linear equation from a nonlinear one? What does it mean for a function to be a "solution," and when can you be sure a solution exists and is unique?

This module provides that vocabulary. It also introduces the direction field — a visual tool that lets you see the qualitative behaviour of solutions without solving a single equation — and the Euler method, the simplest numerical algorithm for generating approximate solutions. Together, these analytical, visual, and computational perspectives form the lens through which the entire course will be viewed.

Example 1: Classification and Verification

Problem. Classify the equation \(y'' + 4y = \sin(x)\) and verify that \(y_p(x) = \frac{1}{3}\sin(x)\) is a particular solution.

Step 1 — Classification

• Type: ODE (one independent variable, \(x\)).
• Order: 2 (highest derivative is \(y''\)).
• Degree: 1 (\(y''\) appears to the first power).
• Linearity: Linear. The equation has the form \(a_2 y'' + a_0 y = g(x)\) with \(a_2=1\), \(a_0=4\), \(g(x)=\sin(x)\), and the unknown \(y\) and its derivatives appear only to the first power with no mutual products.

Step 2 — Verification

Compute derivatives of \(y_p = \frac{1}{3}\sin(x)\):

Substitute into the left side:

This equals the right side, so \(y_p = \frac{1}{3}\sin(x)\) is indeed a particular solution.

Step 3 — General Solution

The associated homogeneous equation \(y'' + 4y = 0\) has characteristic equation \(r^2 + 4 = 0\), giving \(r = \pm 2i\). The homogeneous solution is \(y_h = c_1\cos(2x) + c_2\sin(2x)\). Therefore the general solution is:

Interpretation

The general solution contains two arbitrary constants \(c_1, c_2\) because the equation is second-order. The homogeneous part represents free oscillations at the natural frequency \(\omega_n = 2\), while the particular solution represents forced oscillation at the driving frequency \(\omega = 1\). When the driving frequency equals the natural frequency, resonance occurs — a phenomenon explored in Level 1.

Example 2: Existence, Uniqueness, and Direction Field Analysis

Problem. Consider the IVP \(y' = x^2 + y^2,\quad y(0) = 0\). (a) Does the Picard-Lindelöf theorem guarantee a unique local solution? (b) Describe the direction field qualitatively.

Step 1 — Check Picard-Lindelöf Hypotheses

Here \(f(x,y) = x^2 + y^2\). This function is a polynomial in \(x\) and \(y\), so it is continuous everywhere. Its partial derivative with respect to \(y\) is

which is also continuous everywhere. Both hypotheses of the Picard-Lindelöf theorem are satisfied in any rectangle containing \((0,0)\).

Conclusion: There exists a unique solution in some neighbourhood of \(x=0\).

Step 2 — Direction Field Analysis

The slope at any point \((x,y)\) is \(f(x,y) = x^2 + y^2 \ge 0\). This means:

• All slopes are non-negative. Solution curves never decrease.
• At the origin \((0,0)\), the slope is 0 (horizontal tangent).
• Along the \(y\)-axis (\(x=0\)), the slope is \(y^2\), which is zero at the origin and grows quadratically.
• As \(|x|\) or \(|y|\) increases, slopes become very steep, suggesting rapid growth.

This equation is a Riccati equation and cannot be solved in elementary terms. Its solution is related to Bessel functions and exhibits finite-time blowup: the solution \(y(x)\) goes to \(+\infty\) at a finite value of \(x\). The direction field makes this visible — solution curves curve upward with accelerating steepness.

Step 3 — Euler's Method Approximation

Using step size \(h = 0.1\) starting from \((x_0, y_0) = (0, 0)\):

Even this crude approximation shows the solution beginning to curve upward, consistent with the direction field analysis.

Example 3: Non-Uniqueness When Picard-Lindelöf Fails

Problem. Show that the IVP \(y' = 3y^{2/3},\quad y(0) = 0\) has at least two solutions.

Step 1 — Check the Hypotheses

Here \(f(x,y) = 3y^{2/3}\) is continuous for all \(y\), but

which is undefined (and unbounded) at \(y=0\). The Lipschitz condition fails at the initial point. The theorem does not guarantee uniqueness.

Step 2 — Find Two Solutions

Solution 1: \(y(x) = 0\) for all \(x\). Check: \(y' = 0\) and \(3(0)^{2/3} = 0\). Valid.

Solution 2: Try \(y = (x - c)^3\) for \(x \ge c\), \(y=0\) for \(x < c\). With \(c=0\): let \(y = x^3\). Then \(y' = 3x^2\) and \(3y^{2/3} = 3(x^3)^{2/3} = 3x^2\). Valid for \(x \ge 0\).

Step 3 — Infinitely Many Solutions

In fact, for any \(c \ge 0\), the function

satisfies the IVP. This gives an infinite family of solutions, all passing through the origin with initial value 0.

Interpretation

When the Picard-Lindelöf hypotheses fail, uniqueness is not guaranteed, and physically this means the model is incomplete — additional information is needed to determine which trajectory the system actually follows. This is a critical consideration when modelling real systems.

Lab 1: Direction Field Plotter

This lab plots the direction field for a first-order ODE \(y' = f(x,y)\). Modify the function f to visualise different equations.

import numpy as np
import matplotlib
matplotlib.use('Agg')  # For non-interactive environments
import matplotlib.pyplot as plt

def plot_direction_field(f, x_range, y_range, nx=20, ny=20, title="Direction Field"):
    """
    Plot the direction field of y' = f(x, y).

    Parameters
    ----------
    f : callable
        Function f(x, y) returning the slope y'.
    x_range : tuple
        (x_min, x_max) for the plot domain.
    y_range : tuple
        (y_min, y_max) for the plot range.
    nx, ny : int
        Number of grid points in each direction.
    """
    x = np.linspace(x_range[0], x_range[1], nx)
    y = np.linspace(y_range[0], y_range[1], ny)
    X, Y = np.meshgrid(x, y)

    # Compute slopes
    DY = f(X, Y)
    DX = np.ones_like(DY)

    # Normalise arrow lengths for uniform appearance
    N = np.sqrt(DX**2 + DY**2)
    N[N == 0] = 1  # Avoid division by zero
    DX_norm = DX / N
    DY_norm = DY / N

    fig, ax = plt.subplots(figsize=(8, 6))
    ax.quiver(X, Y, DX_norm, DY_norm, N,
              cmap='coolwarm', angles='xy', scale=30, width=0.003)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(title)
    ax.set_xlim(x_range)
    ax.set_ylim(y_range)
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig('direction_field.png', dpi=150)
    plt.show()
    print("Direction field saved to direction_field.png")

# === Example: y' = x^2 + y^2 ===
f = lambda x, y: x**2 + y**2
plot_direction_field(f, (-2, 2), (-2, 2), title=r"Direction field: $y' = x^2 + y^2$")

# === Try another equation: y' = -y + sin(x) ===
g = lambda x, y: -y + np.sin(x)
plot_direction_field(g, (-4, 4), (-3, 3), title=r"Direction field: $y' = -y + \sin(x)$")

Expected output: Two direction field plots. The first (\(y' = x^2 + y^2\)) shows all arrows pointing upward with slopes growing rapidly away from the origin. The second (\(y' = -y + \sin(x)\)) shows arrows converging toward a sinusoidal attractor, illustrating stable equilibrium behaviour.

Lab 2: Euler's Method with Error Analysis

We implement Euler's method, apply it to the IVP \(y' = -2y,\; y(0)=1\) (whose exact solution is \(y = e^{-2x}\)), and measure how the error depends on the step size \(h\).

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

def euler_method(f, x0, y0, x_end, h):
    """
    Solve y' = f(x, y), y(x0) = y0 using Euler's method.

    Returns arrays of x-values and y-values.
    """
    n_steps = int((x_end - x0) / h)
    x = np.zeros(n_steps + 1)
    y = np.zeros(n_steps + 1)
    x[0], y[0] = x0, y0

    for i in range(n_steps):
        y[i+1] = y[i] + h * f(x[i], y[i])
        x[i+1] = x[i] + h

    return x, y

# Define the ODE and exact solution
f = lambda x, y: -2 * y
exact = lambda x: np.exp(-2 * x)

x0, y0, x_end = 0.0, 1.0, 3.0

# --- Plot solutions for different step sizes ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

x_exact = np.linspace(x0, x_end, 200)
ax1.plot(x_exact, exact(x_exact), 'k-', linewidth=2, label='Exact: $y = e^{-2x}$')

step_sizes = [0.5, 0.2, 0.1, 0.05]
errors_at_end = []

for h in step_sizes:
    x_euler, y_euler = euler_method(f, x0, y0, x_end, h)
    ax1.plot(x_euler, y_euler, 'o--', markersize=3, label=f'Euler h={h}')
    error = abs(y_euler[-1] - exact(x_euler[-1]))
    errors_at_end.append(error)
    print(f"h = {h:.3f}: y({x_end}) = {y_euler[-1]:.6f}, "
          f"exact = {exact(x_end):.6f}, error = {error:.6e}")

ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title("Euler's Method: $y' = -2y$, $y(0) = 1$")
ax1.legend()
ax1.grid(True, alpha=0.3)

# --- Convergence plot ---
ax2.loglog(step_sizes, errors_at_end, 'bo-', linewidth=2, label='Measured error')
ax2.loglog(step_sizes, [s * errors_at_end[0] / step_sizes[0] for s in step_sizes],
           'r--', label='Slope 1 reference')
ax2.set_xlabel('Step size h')
ax2.set_ylabel('|Error at x=3|')
ax2.set_title('Error vs. Step Size (log-log)')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('euler_convergence.png', dpi=150)
plt.show()
print("\nPlot saved to euler_convergence.png")
print("The log-log slope is approximately 1, confirming first-order convergence.")

Expected output: The left panel shows the exact exponential decay alongside Euler approximations that improve as \(h\) decreases. The right panel shows that error vs. step size follows a straight line with slope 1 on a log-log plot, confirming that Euler's method is first-order accurate.

Exercise A1. Classify each equation by order, degree, type (ODE/PDE), and linearity:

1. \(y''' - 2y' + y = e^x\)
2. \((y')^2 + y = x\)
3. \(\dfrac{\partial^2 u}{\partial x^2} + \dfrac{\partial^2 u}{\partial y^2} = 0\)
4. \(y\,y'' + (y')^2 = 0\)

Exercise A2. Verify that \(y(x) = c_1 e^x + c_2 e^{-x}\) is the general solution of \(y'' - y = 0\) for arbitrary constants \(c_1, c_2\). Then find the particular solution satisfying \(y(0) = 2\), \(y'(0) = 0\).

Exercise A3. For the IVP \(y' = y^2,\; y(0) = 1\), find the exact solution by separation of variables. At what value of \(x\) does the solution blow up (become infinite)?

Exercise A4. Does the Picard-Lindelöf theorem guarantee a unique solution for the IVP \(y' = \frac{x}{y},\; y(0) = 0\)? Justify your answer carefully.

Exercise C1. Implement a direction-field plotter for \(y' = y - x^2 + 1\) on the domain \([-1, 4] \times [-1, 4]\). Overlay three Euler-method solution curves starting from \(y(0) = 0\), \(y(0) = 0.5\), and \(y(0) = 2\), all with step size \(h = 0.05\). Do the trajectories converge or diverge?

Exercise C2. For the IVP \(y' = -2y + e^{-x},\; y(0) = 1\), the exact solution is \(y = e^{-x} + 0\cdot e^{-2x} = e^{-x}\) (verify this). Run Euler's method with \(h = 0.1, 0.01, 0.001\) on \([0, 5]\). For each \(h\), compute the maximum absolute error over the interval and verify that the error scales linearly with \(h\).

import numpy as np

def euler(f, x0, y0, x_end, h):
    n = int((x_end - x0) / h)
    x = np.linspace(x0, x_end, n+1)
    y = np.zeros(n+1)
    y[0] = y0
    for i in range(n):
        y[i+1] = y[i] + h * f(x[i], y[i])
    return x, y

f = lambda x, y: -2*y + np.exp(-x)
exact = lambda x: np.exp(-x)

for h in [0.1, 0.01, 0.001]:
    x, y = euler(f, 0, 1, 5, h)
    max_err = np.max(np.abs(y - exact(x)))
    print(f"h = {h:.4f}: max error = {max_err:.6e}")
# Expected: errors decrease by factor ~10 each time (first-order method).

Exercise C3. Write a Python function that, given \(f(x,y)\) and a point \((x_0, y_0)\), checks the Picard-Lindelöf conditions numerically. Specifically, your function should: (a) evaluate \(f\) at and near the point to check continuity (using finite differences), (b) approximate \(\partial f / \partial y\) using a central difference, and (c) check whether the partial derivative appears bounded near the point. Test on \(f(x,y) = 3y^{2/3}\) at \((0,0)\) and on \(f(x,y) = x^2 + y^2\) at \((0,0)\).

Exercise C4. Implement the improved Euler method (Heun's method):

Apply it to \(y' = -2y,\; y(0) = 1\) on \([0, 3]\). Compare the error with standard Euler for step sizes \(h = 0.5, 0.2, 0.1, 0.05\). Verify that the improved Euler method is second-order (error scales as \(h^2\)).

import numpy as np

def improved_euler(f, x0, y0, x_end, h):
    n = int((x_end - x0) / h)
    x = np.linspace(x0, x_end, n+1)
    y = np.zeros(n+1)
    y[0] = y0
    for i in range(n):
        k1 = h * f(x[i], y[i])
        k2 = h * f(x[i] + h, y[i] + k1)
        y[i+1] = y[i] + 0.5 * (k1 + k2)
    return x, y

f = lambda x, y: -2 * y
exact_end = np.exp(-6)  # y(3)

for h in [0.5, 0.2, 0.1, 0.05]:
    _, ye = improved_euler(f, 0, 1, 3, h)
    err = abs(ye[-1] - exact_end)
    print(f"h = {h:.2f}: error = {err:.6e}")
# Ratio of errors when h halves should be ~4 (second-order).

Exercise G1. (Policy comparison.) Consider the ODE \(y' = \sin(x) - y\). Treat the direction field as a policy map. Use Python to generate the direction field and overlay two solution trajectories: one starting at \((0, -2)\) and one at \((0, 3)\). Both converge to the same long-term behaviour. In RL terms, this corresponds to a policy with a single attracting fixed point (or cycle). Determine the asymptotic (long-time) behaviour analytically by guessing a particular solution of the form \(y_p = A\sin(x) + B\cos(x)\) and finding \(A\) and \(B\). Interpret: why does every initial condition lead to the same eventual trajectory?

Exercise G2. (Exploration and finite-time blowup.) An agent following the policy \(y' = 1 + y^2\) (with \(y(0)=0\)) will reach \(y = +\infty\) in finite time. (a) Solve the ODE exactly. (b) What is the blowup time? (c) Run Euler's method with \(h=0.01\). At what step does the numerical solution first exceed \(10^6\)? (d) In RL language, this is a policy that leads the agent "off a cliff." Design a modified policy \(\tilde{f}(x,y) = \min(1+y^2,\; M)\) with a reward-clipping parameter \(M\). How does the clipped trajectory differ? Plot both.

1. The equation \(y''' + xy' - y^2 = \ln(x)\) is:
   (a) 3rd order, linear   (b) 3rd order, nonlinear   (c) 2nd order, nonlinear   (d) 3rd order, degree 2
2. Which of the following is a PDE?
   (a) \(y'' + y = 0\)   (b) \(\frac{du}{dt} = -ku\)   (c) \(u_{xx} + u_{yy} = 0\)   (d) \(\frac{dy}{dx} = xy\)
3. The general solution of a 3rd-order ODE contains how many arbitrary constants?
   (a) 1   (b) 2   (c) 3   (d) It depends on whether the equation is linear
4. To verify that \(y = e^{3x}\) is a solution of \(y' - 3y = 0\), you should:
   (a) Integrate \(y' - 3y\)   (b) Compute \(y' = 3e^{3x}\) and check that \(3e^{3x} - 3e^{3x} = 0\)   (c) Find the general solution first   (d) Plot the direction field
5. The IVP \(y' = y^{1/3},\; y(0) = 0\) fails which hypothesis of Picard-Lindelöf?
   (a) \(f\) is not continuous   (b) \(\partial f/\partial y\) is not continuous at the initial point   (c) Both hypotheses fail   (d) Neither hypothesis fails
6. A direction field for \(y' = -y\) shows:
   (a) Horizontal lines everywhere   (b) Slopes that are negative when \(y > 0\) and positive when \(y < 0\)   (c) Vertical lines everywhere   (d) Slopes independent of \(y\)
7. Euler's method with step size \(h\) has global error of order:
   (a) \(h^2\)   (b) \(h\)   (c) \(h^{1/2}\)   (d) \(h^3\)
8. A singular solution of a DE is one that:
   (a) Blows up in finite time   (b) Cannot be obtained from the general solution for any value of the constants   (c) Satisfies the homogeneous equation   (d) Has no derivatives
9. The equation \(y' = x + y\) with \(y(0) = 1\) satisfies the Picard-Lindelöf conditions because:
   (a) The equation is linear   (b) \(f(x,y) = x+y\) and \(\partial f/\partial y = 1\) are both continuous everywhere   (c) The solution can be found by separation of variables   (d) The direction field has no vertical tangents
10. If \(y_1\) and \(y_2\) are solutions of a nonlinear ODE, then \(y_1 + y_2\) is:
    (a) Always a solution   (b) Never a solution   (c) A solution only if \(y_1 = -y_2\)   (d) Not guaranteed to be a solution

Answer Key

1: (b) — \(y^2\) makes it nonlinear   2: (c)   3: (c)   4: (b)   5: (b) — \(\partial f/\partial y = \frac{1}{3}y^{-2/3}\) is undefined at \(y=0\)   6: (b)   7: (b)   8: (b)   9: (b)   10: (d)

Project Description

Build a Python script (or Jupyter notebook) that serves as an interactive exploration tool for first-order ODEs. Your tool should:

1. Accept a user-specified function \(f(x,y)\) (as a Python lambda or string expression) and domain bounds.
2. Plot the direction field.
3. Overlay Euler-method solution curves for user-specified initial conditions (at least 5 different initial conditions on one plot).
4. Compute and display the Euler-method approximation error at \(x=x_{\text{end}}\) for a user-specified exact solution (if known).
5. Automatically check the Picard-Lindelöf conditions at each initial point: evaluate \(\partial f/\partial y\) numerically and warn if it appears unbounded.

Test your tool on the following three ODEs:

• \(y' = -y + \sin(x)\) with initial conditions \(y(0) \in \{-2, -1, 0, 1, 2\}\).
• \(y' = y^2 - x\) with initial conditions \(y(0) \in \{-1, -0.5, 0, 0.5, 1\}\).
• \(y' = \sqrt{|y|}\) with initial condition \(y(0) = 0\) (the Picard-Lindelöf warning should trigger).

Rubric (30 points)