Ordinary Differential Equations

Module 2 — Second-Order Linear ODEs (Constant Coefficients)

def characteristic_roots(a_coef, b_coef, c_coef):
    """
    Characteristic equation for a*y'' + b*y' + c*y = 0:
    ar^2 + br + c = 0  =>  r = (-b +/- sqrt(b^2 - 4ac)) / (2a)

    Discriminant determines solution form:
      D > 0: Two real distinct roots r1, r2
             y_h = C1*exp(r1*t) + C2*exp(r2*t)
      D = 0: Repeated root r = -b/(2a)
             y_h = (C1 + C2*t)*exp(r*t)
      D < 0: Complex conjugate roots alpha +/- i*beta
             y_h = exp(alpha*t) * (C1*cos(beta*t) + C2*sin(beta*t))
    """
    D = b_coef**2 - 4.0 * a_coef * c_coef

    if D > 1e-12:  # Two real distinct roots
        r1 = (-b_coef + math.sqrt(D)) / (2.0 * a_coef)
        r2 = (-b_coef - math.sqrt(D)) / (2.0 * a_coef)
        return {
            "discriminant": round(D, 8),
            "root_type": "real distinct",
            "r1": round(r1, 8), "r2": round(r2, 8),
            "y_h": "C1*exp({:.4f}*t) + C2*exp({:.4f}*t)".format(r1, r2),
            "behavior": "overdamped" if r1 < 0 and r2 < 0 else "unstable" if r1 > 0 or r2 > 0 else "mixed"
        }
    elif abs(D) <= 1e-12:  # Repeated root
        r = -b_coef / (2.0 * a_coef)
        return {
            "discriminant": 0.0,
            "root_type": "repeated real",
            "r": round(r, 8),
            "y_h": "(C1 + C2*t)*exp({:.4f}*t)".format(r),
            "behavior": "critically damped" if r < 0 else "linear growth"
        }
    else:  # Complex roots
        alpha = -b_coef / (2.0 * a_coef)
        beta = math.sqrt(-D) / (2.0 * a_coef)
        return {
            "discriminant": round(D, 8),
            "root_type": "complex conjugate",
            "alpha": round(alpha, 8), "beta": round(beta, 8),
            "y_h": "exp({:.4f}*t) * (C1*cos({:.4f}*t) + C2*sin({:.4f}*t))".format(alpha, beta, beta),
            "omega_n": round(beta, 8),
            "behavior": "underdamped oscillation" if alpha < 0 else ("undamped" if abs(alpha) < 1e-10 else "oscillatory growth")
        }

def undetermined_coefficients_guess(forcing_type, coefficients=None):
    """
    Undetermined coefficients method: guess form for y_p based on right-hand side.
    forcing_type: 'polynomial', 'exponential', 'sinusoidal', 'poly_exp', 'poly_sin'
    Returns: guess form for particular solution
    """
    table = {
        "polynomial_n": "y_p = A_n*t^n + A_{n-1}*t^(n-1) + ... + A_1*t + A_0",
        "exponential":  "y_p = A*exp(alpha*t)  [if alpha is NOT a root of char. eq.]",
        "exponential_resonance": "y_p = A*t*exp(alpha*t)  [if alpha IS a root]",
        "sinusoidal":   "y_p = A*cos(omega*t) + B*sin(omega*t)  [both terms always]",
        "sin_resonance": "y_p = t*(A*cos(omega*t) + B*sin(omega*t))  [if omega matches char. freq.]",
        "poly_exp":     "y_p = (A_n*t^n + ... + A_0) * exp(alpha*t)",
        "poly_sin":     "y_p = (A_n*t^n + ... + A_0)*cos(omega*t) + (B_n*t^n + ...)*sin(omega*t)"
    }
    return table.get(forcing_type, {k: v for k, v in table.items()})

def second_order_constant_coeff_ivp(a, b, c, f_type, f_params,
                                     y0, dy0, t_values):
    """
    Solve a*y'' + b*y' + c*y = 0 numerically using RK4 for general IVP.
    Converts to system: y1 = y, y2 = y'
      y1' = y2
      y2' = (f(t) - b*y2 - c*y1) / a
    f_type: 'zero' (homogeneous), 'const', 'sin', 'exp'
    f_params: dict with 'F' (amplitude), 'omega', 'alpha' as needed
    Returns: list of (t, y, y') values
    """
    def forcing(t):
        if f_type == 'zero':
            return 0.0
        elif f_type == 'const':
            return f_params.get('F', 0.0)
        elif f_type == 'sin':
            return f_params.get('F', 1.0) * math.sin(f_params.get('omega', 1.0) * t)
        elif f_type == 'cos':
            return f_params.get('F', 1.0) * math.cos(f_params.get('omega', 1.0) * t)
        elif f_type == 'exp':
            return f_params.get('F', 1.0) * math.exp(f_params.get('alpha', -1.0) * t)
        return 0.0

    def rhs(t, Y):
        y1, y2 = Y
        return [y2, (forcing(t) - b*y2 - c*y1) / a]

    result = [(t_values[0], y0, dy0)]
    Y = [y0, dy0]

    for i in range(1, len(t_values)):
        h = t_values[i] - t_values[i-1]
        t = t_values[i-1]
        k1 = rhs(t, Y)
        k2 = rhs(t + h/2, [Y[j] + h/2*k1[j] for j in range(2)])
        k3 = rhs(t + h/2, [Y[j] + h/2*k2[j] for j in range(2)])
        k4 = rhs(t + h, [Y[j] + h*k3[j] for j in range(2)])
        Y = [Y[j] + h/6*(k1[j]+2*k2[j]+2*k3[j]+k4[j]) for j in range(2)]
        result.append((round(t_values[i], 8), round(Y[0], 8), round(Y[1], 8)))

    return result

Module 3 — Laplace Transforms

def laplace_table():
    """
    Standard Laplace transform pairs L{f(t)} = F(s).
    All pairs valid for t >= 0.
    """
    return {
        "1":           "1/s",
        "t":           "1/s^2",
        "t^n":         "n! / s^(n+1)",
        "e^(at)":      "1/(s-a)   [s > a]",
        "t*e^(at)":    "1/(s-a)^2",
        "t^n*e^(at)":  "n!/(s-a)^(n+1)",
        "sin(bt)":     "b/(s^2+b^2)",
        "cos(bt)":     "s/(s^2+b^2)",
        "e^(at)*sin(bt)": "b/((s-a)^2+b^2)",
        "e^(at)*cos(bt)": "(s-a)/((s-a)^2+b^2)",
        "t*sin(bt)":   "2bs/(s^2+b^2)^2",
        "t*cos(bt)":   "(s^2-b^2)/(s^2+b^2)^2",
        "u(t-c)":      "e^(-cs)/s   [unit step at t=c]",
        "delta(t-c)":  "e^(-cs)   [Dirac delta at t=c]",
        "f'(t)":       "s*F(s) - f(0)   [derivative property]",
        "f''(t)":      "s^2*F(s) - s*f(0) - f'(0)   [2nd derivative]",
    }

def laplace_derivative_property(s, F_s_expr, y0, dy0=None, order=1):
    """
    Laplace transform of derivatives:
    L{y'} = s*Y(s) - y(0)
    L{y''} = s^2*Y(s) - s*y(0) - y'(0)
    Returns: transformed expression (symbolic guidance)
    """
    if order == 1:
        return {"L_y_prime": f"s*Y(s) - {y0}"}
    elif order == 2:
        return {"L_y_double_prime": f"s^2*Y(s) - {y0}*s - {dy0}"}
    return None

def inverse_laplace_partial_fractions(numerator_str, denominator_str):
    """
    Guide for partial fraction decomposition before inverse Laplace.
    Denominator root forms and their inverse transforms:
    """
    return {
        "method": "Partial Fraction Decomposition",
        "forms": {
            "A/(s-a)":             "A*e^(a*t)",
            "A/(s-a)^2":           "A*t*e^(a*t)",
            "A/(s-a)^n":           "A/(n-1)! * t^(n-1)*e^(a*t)",
            "(As+B)/(s^2+b^2)":    "A*cos(bt) + (B/b)*sin(bt)",
            "(As+B)/(s^2+2as+a^2+b^2)": "e^(-at)*[A*cos(bt) + C*sin(bt)]",
        },
        "steps": [
            "1. Factor denominator completely into linear and irreducible quadratic factors",
            "2. Write partial fraction expansion with unknown numerators",
            "3. Multiply through and match coefficients (or use Heaviside cover-up for simple poles)",
            "4. Look up each term in Laplace table",
        ],
        "tip": "Complete the square for quadratics: s^2+bs+c = (s+b/2)^2 + (c - b^2/4)"
    }

Module 4 — Numerical Methods (Euler and RK4)

def euler_method(f, t0, y0, h, n_steps):
    """
    Euler's method: y_{n+1} = y_n + h * f(t_n, y_n)
    f:       callable f(t, y) -> dy/dt
    t0, y0:  initial conditions
    h:       step size
    n_steps: number of steps
    Returns: list of (t, y) pairs
    First-order, O(h) global error — use RK4 for better accuracy.
    """
    result = [(t0, y0)]
    t, y = t0, y0
    for _ in range(n_steps):
        y = y + h * f(t, y)
        t = t + h
        result.append((round(t, 10), round(y, 10)))
    return result

def rk4_method(f, t0, y0, h, n_steps):
    """
    4th-order Runge-Kutta: O(h^4) global error, much more accurate than Euler.
    k1 = h*f(t, y)
    k2 = h*f(t+h/2, y+k1/2)
    k3 = h*f(t+h/2, y+k2/2)
    k4 = h*f(t+h, y+k3)
    y_{n+1} = y_n + (k1 + 2k2 + 2k3 + k4)/6
    f:       callable f(t, y) -> dy/dt
    Returns: list of (t, y) pairs
    """
    result = [(t0, y0)]
    t, y = t0, y0
    for _ in range(n_steps):
        k1 = h * f(t, y)
        k2 = h * f(t + h/2.0, y + k1/2.0)
        k3 = h * f(t + h/2.0, y + k2/2.0)
        k4 = h * f(t + h, y + k3)
        y = y + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
        t = t + h
        result.append((round(t, 10), round(y, 10)))
    return result

def rk4_system(f_list, t0, y0_list, h, n_steps):
    """
    RK4 for a system of ODEs: dy_i/dt = f_i(t, y_1, ..., y_n)
    f_list:   list of callables [f1, f2, ..., fn] where f_i(t, y_list) -> dy_i/dt
    y0_list:  list of initial conditions [y1(0), y2(0), ..., yn(0)]
    Returns: list of (t, y_1, y_2, ..., y_n) tuples
    """
    n = len(f_list)
    result = [(t0,) + tuple(y0_list)]
    t = t0
    Y = list(y0_list)

    for _ in range(n_steps):
        k1 = [h * f_list[i](t, Y) for i in range(n)]
        Y2 = [Y[i] + k1[i]/2.0 for i in range(n)]
        k2 = [h * f_list[i](t + h/2.0, Y2) for i in range(n)]
        Y3 = [Y[i] + k2[i]/2.0 for i in range(n)]
        k3 = [h * f_list[i](t + h/2.0, Y3) for i in range(n)]
        Y4 = [Y[i] + k3[i] for i in range(n)]
        k4 = [h * f_list[i](t + h, Y4) for i in range(n)]
        Y = [Y[i] + (k1[i] + 2.0*k2[i] + 2.0*k3[i] + k4[i])/6.0 for i in range(n)]
        t = t + h
        result.append((round(t, 8),) + tuple(round(y, 8) for y in Y))

    return result

Module 5 — Systems of ODEs and Phase Plane

def eigenvalue_2x2(a11, a12, a21, a22):
    """
    Eigenvalues of 2x2 matrix [[a11, a12], [a21, a22]]:
    lambda^2 - trace*lambda + det = 0
    trace = a11 + a22,  det = a11*a22 - a12*a21
    Returns: eigenvalues, eigenvector hints, phase plane type
    """
    trace = a11 + a22
    det = a11*a22 - a12*a21
    discriminant = trace**2 - 4.0*det

    if discriminant > 1e-10:
        l1 = (trace + math.sqrt(discriminant)) / 2.0
        l2 = (trace - math.sqrt(discriminant)) / 2.0
        eigenvalues = [round(l1, 8), round(l2, 8)]
    elif abs(discriminant) <= 1e-10:
        l = trace / 2.0
        eigenvalues = [round(l, 8), round(l, 8)]
    else:
        alpha = trace / 2.0
        beta = math.sqrt(-discriminant) / 2.0
        eigenvalues = [f"{round(alpha,5)} + {round(beta,5)}i",
                       f"{round(alpha,5)} - {round(beta,5)}i"]
        alpha_val = alpha
        beta_val = beta

    result = {"trace": round(trace, 6), "det": round(det, 6),
              "eigenvalues": eigenvalues}
    result.update(phase_plane_type(trace, det, discriminant))
    return result

def phase_plane_type(trace, det, discriminant=None):
    """
    Classify the equilibrium point for a linear 2D system x' = Ax.
    Based on trace and determinant of A:
    """
    if det < 0:
        return {"equilibrium_type": "saddle point", "stability": "unstable",
                "behavior": "Solutions diverge along one eigenvector, converge along other"}
    elif det == 0:
        return {"equilibrium_type": "degenerate (line of equilibria)"}
    elif discriminant is not None and discriminant > 0:
        if trace < 0:
            return {"equilibrium_type": "stable node (nodal sink)",
                    "stability": "asymptotically stable",
                    "behavior": "All trajectories approach origin as t->inf"}
        elif trace > 0:
            return {"equilibrium_type": "unstable node (nodal source)",
                    "stability": "unstable",
                    "behavior": "All trajectories diverge from origin"}
        else:
            return {"equilibrium_type": "center (star node)", "stability": "neutral"}
    elif discriminant is not None and discriminant < 0:
        if trace < 0:
            return {"equilibrium_type": "stable spiral (spiral sink)",
                    "stability": "asymptotically stable",
                    "behavior": "Spirals inward to origin; damped oscillation"}
        elif trace > 0:
            return {"equilibrium_type": "unstable spiral (spiral source)",
                    "stability": "unstable"}
        else:
            return {"equilibrium_type": "center",
                    "stability": "neutrally stable (periodic orbits)"}
    return {"equilibrium_type": "complex or degenerate; check manually"}

Learning Resources

Textbooks

Free Online Resources (All Verified High-Quality)

WVU Courses

• MATH 261 Applied Differential Equations — first/second order ODEs, systems, Laplace
• MATH 441 Applied Math / PDEs — partial differential equations (graduate level)
• Related: ChBE 311/321 use ODE models; PNGE 321 uses material balance ODEs

Quick Method Reference

Output Format

## ODE Solution — [Type]: [Problem Description]

### ODE Classification
Type: [1st order linear / 2nd order constant-coeff / system / etc.]
Method: [separation / integrating factor / char. eq. / Laplace / RK4]

### Solution Steps
Step 1 — [Identify form]:
[write ODE in standard form]

Step 2 — [Apply method]:
[show calculation]

Step 3 — [General solution]:
y_h = [homogeneous solution]
y_p = [particular solution, if applicable]
y = y_h + y_p = [general solution]

Step 4 — [Apply IVP if given]:
y(0) = [value] => C1 = [value]
y'(0) = [value] => C2 = [value]

### Final Solution
y(t) = [explicit formula]

### Numerical Check (RK4, if analytical verification needed)
| t | y (analytical) | y (RK4) | error |
|---|---------------|---------|-------|

**Behavior:** [stable/unstable/oscillatory/growing — describe physical interpretation]

Error Handling

Caveats

• Numerical stability (stiff ODEs): If the ODE has eigenvalues with vastly different magnitudes (e.g., fast chemical reactions with slow temperature change), explicit methods (Euler, RK4) require very small h. Use implicit methods or scipy.integrate.solve_ivp with method='Radau' or 'BDF' in practice.
• Phase plane analysis is for LINEAR systems only (or linearized near equilibrium). Nonlinear systems can have limit cycles and other features not predicted by linearization.
• Undetermined coefficients vs. Variation of Parameters: UndCoeff works only when g(t) is a polynomial, exponential, sine/cosine, or product of these. Variation of Parameters (Fogler) is more general but more tedious.
• Laplace transform requires zero initial conditions for standard tables. If y(0) != 0 or y'(0) != 0, include the extra s*y(0) and y'(0) terms from the derivative property.
• For PDEs (heat equation, wave equation, Laplace equation) — which appear in ChBE and PNGE transport problems — MATH 261 introduces the separation of variables technique; MATH 441 covers this in depth.