differential_eq.tex

\documentclass[a4paper, 11pt]{report}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{xcolor}
\usepackage[makestderr]{pythontex}

\definecolor{stable}{RGB}{46, 204, 113}
\definecolor{unstable}{RGB}{231, 76, 60}
\definecolor{saddle}{RGB}{241, 196, 15}
\definecolor{limit}{RGB}{155, 89, 182}

\title{Ordinary Differential Equations:\\
Phase Portraits, Stability Analysis, and Limit Cycles}
\author{Department of Applied Mathematics\\Technical Report AM-2024-002}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This report provides a comprehensive computational analysis of ordinary differential equations (ODEs). We examine first and second-order ODEs, construct phase portraits for autonomous systems, perform stability analysis of equilibrium points, and investigate limit cycles in nonlinear oscillators. The van der Pol oscillator is analyzed in detail to demonstrate relaxation oscillations and Hopf bifurcations.
\end{abstract}

\tableofcontents

\chapter{Introduction}

Ordinary differential equations describe the evolution of systems in terms of derivatives with respect to a single variable. This report focuses on qualitative analysis through phase portraits and stability theory.

\section{Classification of ODEs}
\begin{itemize}
    \item \textbf{First-order}: $\frac{dy}{dt} = f(t, y)$
    \item \textbf{Second-order}: $\frac{d^2y}{dt^2} = f(t, y, \frac{dy}{dt})$
    \item \textbf{Linear}: Coefficients are functions of $t$ only
    \item \textbf{Autonomous}: No explicit time dependence
\end{itemize}

\chapter{First-Order ODEs}

\section{Analytical Solutions}
Consider the first-order linear ODE:
\begin{equation}
\frac{dy}{dt} + p(t)y = q(t)
\end{equation}

Solution via integrating factor $\mu(t) = e^{\int p(t)dt}$:
\begin{equation}
y(t) = \frac{1}{\mu(t)}\left[\int \mu(t)q(t)dt + C\right]
\end{equation}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, odeint
from scipy.optimize import fsolve
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

# First-order examples
def exponential_decay(t, y, k=1.0):
    return -k * y

def logistic_growth(t, y, r=1.0, K=10.0):
    return r * y * (1 - y/K)

# Solve multiple first-order ODEs
t_span = (0, 10)
t_eval = np.linspace(0, 10, 200)

# Exponential decay with different initial conditions
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

ax = axes[0, 0]
for y0 in [1, 2, 3, 4, 5]:
    sol = solve_ivp(exponential_decay, t_span, [y0], t_eval=t_eval, args=(0.5,))
    ax.plot(sol.t, sol.y[0], label=f'$y_0 = {y0}$')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$y(t)$')
ax.set_title('Exponential Decay: $dy/dt = -ky$')
ax.legend()
ax.grid(True, alpha=0.3)

# Logistic growth
ax = axes[0, 1]
for y0 in [0.5, 2, 5, 8, 12, 15]:
    sol = solve_ivp(logistic_growth, t_span, [y0], t_eval=t_eval)
    ax.plot(sol.t, sol.y[0], label=f'$y_0 = {y0}$')
ax.axhline(10, color='red', linestyle='--', alpha=0.5, label='$K = 10$')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$y(t)$')
ax.set_title('Logistic Growth: $dy/dt = ry(1-y/K)$')
ax.legend(loc='right')
ax.grid(True, alpha=0.3)

# Direction field for dy/dt = y - y^2
ax = axes[1, 0]
y_range = np.linspace(-0.5, 1.5, 20)
t_range = np.linspace(0, 5, 20)
T, Y = np.meshgrid(t_range, y_range)
dY = Y - Y**2
dT = np.ones_like(dY)
norm = np.sqrt(dT**2 + dY**2)
ax.quiver(T, Y, dT/norm, dY/norm, alpha=0.5)

# Plot solutions
for y0 in [0.1, 0.5, 0.9, 1.1, 1.5]:
    sol = solve_ivp(lambda t, y: y - y**2, (0, 5), [y0], t_eval=np.linspace(0, 5, 100))
    ax.plot(sol.t, sol.y[0], linewidth=2)
ax.set_xlabel('Time $t$')
ax.set_ylabel('$y(t)$')
ax.set_title('Direction Field: $dy/dt = y - y^2$')
ax.set_xlim(0, 5)
ax.set_ylim(-0.5, 1.5)
ax.grid(True, alpha=0.3)

# Bernoulli equation
ax = axes[1, 1]
# dy/dt = y - y^3 (bistable)
for y0 in np.linspace(-1.5, 1.5, 15):
    sol = solve_ivp(lambda t, y: y - y**3, (0, 5), [y0], t_eval=np.linspace(0, 5, 100))
    ax.plot(sol.t, sol.y[0], 'b-', linewidth=0.8, alpha=0.6)
ax.axhline(1, color='green', linestyle='--', label='Stable: $y=1$')
ax.axhline(-1, color='green', linestyle='--', label='Stable: $y=-1$')
ax.axhline(0, color='red', linestyle='--', label='Unstable: $y=0$')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$y(t)$')
ax.set_title('Bistable System: $dy/dt = y - y^3$')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('first_order_odes.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{first_order_odes.pdf}
\caption{First-order ODE solutions: (a) exponential decay, (b) logistic growth, (c) direction field, (d) bistable system.}
\end{figure}

\chapter{Second-Order Linear ODEs}

\section{Harmonic Oscillator}
The damped harmonic oscillator:
\begin{equation}
m\frac{d^2x}{dt^2} + c\frac{dx}{dt} + kx = 0
\end{equation}

Characteristic equation: $mr^2 + cr + k = 0$

\begin{pycode}
# Second-order ODE as system of first-order
def harmonic_oscillator(t, state, m=1.0, c=0.0, k=1.0):
    x, v = state
    return [v, -(c*v + k*x)/m]

fig, axes = plt.subplots(2, 3, figsize=(14, 8))

# Undamped
t_span = (0, 20)
t_eval = np.linspace(0, 20, 500)

ax = axes[0, 0]
sol = solve_ivp(harmonic_oscillator, t_span, [1, 0], t_eval=t_eval, args=(1, 0, 1))
ax.plot(sol.t, sol.y[0], 'b-')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$x(t)$')
ax.set_title('Undamped ($\\zeta = 0$)')
ax.grid(True, alpha=0.3)

# Underdamped
ax = axes[0, 1]
sol = solve_ivp(harmonic_oscillator, t_span, [1, 0], t_eval=t_eval, args=(1, 0.3, 1))
ax.plot(sol.t, sol.y[0], 'b-')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$x(t)$')
ax.set_title('Underdamped ($\\zeta = 0.15$)')
ax.grid(True, alpha=0.3)

# Critically damped
ax = axes[0, 2]
sol = solve_ivp(harmonic_oscillator, t_span, [1, 0], t_eval=t_eval, args=(1, 2, 1))
ax.plot(sol.t, sol.y[0], 'b-')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$x(t)$')
ax.set_title('Critically Damped ($\\zeta = 1$)')
ax.grid(True, alpha=0.3)

# Overdamped
ax = axes[1, 0]
sol = solve_ivp(harmonic_oscillator, t_span, [1, 0], t_eval=t_eval, args=(1, 4, 1))
ax.plot(sol.t, sol.y[0], 'b-')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$x(t)$')
ax.set_title('Overdamped ($\\zeta = 2$)')
ax.grid(True, alpha=0.3)

# Phase portraits
ax = axes[1, 1]
for zeta in [0, 0.15, 0.5, 1.0]:
    c = 2 * zeta
    sol = solve_ivp(harmonic_oscillator, (0, 30), [1, 0], t_eval=np.linspace(0, 30, 1000), args=(1, c, 1))
    ax.plot(sol.y[0], sol.y[1], label=f'$\\zeta = {zeta}$')
ax.set_xlabel('$x$')
ax.set_ylabel('$\\dot{x}$')
ax.set_title('Phase Portraits')
ax.legend()
ax.grid(True, alpha=0.3)

# Energy decay
ax = axes[1, 2]
for zeta in [0.1, 0.3, 0.5]:
    c = 2 * zeta
    sol = solve_ivp(harmonic_oscillator, t_span, [1, 0], t_eval=t_eval, args=(1, c, 1))
    energy = 0.5 * sol.y[0]**2 + 0.5 * sol.y[1]**2
    ax.plot(sol.t, energy, label=f'$\\zeta = {zeta}$')
ax.set_xlabel('Time $t$')
ax.set_ylabel('Energy $E$')
ax.set_title('Energy Decay')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_yscale('log')

plt.tight_layout()
plt.savefig('second_order_odes.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{second_order_odes.pdf}
\caption{Second-order ODE: damped harmonic oscillator with various damping ratios.}
\end{figure}

\chapter{Phase Portrait Analysis}

\section{2D Autonomous Systems}
Consider the system:
\begin{align}
\frac{dx}{dt} &= f(x, y) \\
\frac{dy}{dt} &= g(x, y)
\end{align}

Equilibrium points satisfy $f(x^*, y^*) = g(x^*, y^*) = 0$.

\section{Linear Stability Analysis}
Near an equilibrium, linearize using the Jacobian:
\begin{equation}
J = \begin{pmatrix}
\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\
\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y}
\end{pmatrix}
\end{equation}

Eigenvalues $\lambda$ determine stability.

\begin{pycode}
# 2D systems with different equilibrium types
def saddle_system(t, state):
    x, y = state
    return [x, -y]

def stable_node(t, state):
    x, y = state
    return [-x, -2*y]

def unstable_spiral(t, state):
    x, y = state
    return [0.1*x - y, x + 0.1*y]

def center(t, state):
    x, y = state
    return [-y, x]

fig, axes = plt.subplots(2, 2, figsize=(12, 12))
systems = [
    (saddle_system, 'Saddle Point'),
    (stable_node, 'Stable Node'),
    (unstable_spiral, 'Unstable Spiral'),
    (center, 'Center')
]

for ax, (system, title) in zip(axes.flatten(), systems):
    # Vector field
    x_range = np.linspace(-3, 3, 20)
    y_range = np.linspace(-3, 3, 20)
    X, Y = np.meshgrid(x_range, y_range)

    state = np.array([X, Y])
    derivatives = system(0, state)
    U, V = derivatives[0], derivatives[1]

    norm = np.sqrt(U**2 + V**2)
    norm[norm == 0] = 1
    ax.quiver(X, Y, U/norm, V/norm, alpha=0.4)

    # Trajectories
    angles = np.linspace(0, 2*np.pi, 8, endpoint=False)
    for angle in angles:
        x0 = 2.5 * np.cos(angle)
        y0 = 2.5 * np.sin(angle)
        try:
            sol = solve_ivp(system, (0, 10), [x0, y0], t_eval=np.linspace(0, 10, 500), max_step=0.1)
            ax.plot(sol.y[0], sol.y[1], 'b-', linewidth=1)
        except:
            pass

    ax.plot(0, 0, 'ro', markersize=8)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_title(title)
    ax.set_xlim(-3, 3)
    ax.set_ylim(-3, 3)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('phase_portraits.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{phase_portraits.pdf}
\caption{Phase portraits showing different equilibrium types based on eigenvalue classification.}
\end{figure}

\chapter{Limit Cycles}

\section{Van der Pol Oscillator}
The van der Pol oscillator exhibits self-sustained oscillations:
\begin{equation}
\frac{d^2x}{dt^2} - \mu(1 - x^2)\frac{dx}{dt} + x = 0
\end{equation}

\begin{pycode}
def van_der_pol(t, state, mu=1.0):
    x, v = state
    return [v, mu*(1 - x**2)*v - x]

fig, axes = plt.subplots(2, 3, figsize=(14, 8))

# Time series for different mu
mu_values = [0.1, 1.0, 5.0]
t_span = (0, 50)
t_eval = np.linspace(0, 50, 2000)

for i, mu in enumerate(mu_values):
    ax = axes[0, i]
    sol = solve_ivp(van_der_pol, t_span, [0.1, 0], t_eval=t_eval, args=(mu,))
    ax.plot(sol.t, sol.y[0], 'b-', linewidth=0.8)
    ax.set_xlabel('Time $t$')
    ax.set_ylabel('$x(t)$')
    ax.set_title(f'Van der Pol: $\\mu = {mu}$')
    ax.grid(True, alpha=0.3)

# Phase portraits with limit cycle
for i, mu in enumerate(mu_values):
    ax = axes[1, i]

    # Multiple initial conditions
    for x0, v0 in [(0.1, 0), (4, 0), (0, 4), (-3, -3)]:
        sol = solve_ivp(van_der_pol, (0, 100), [x0, v0], t_eval=np.linspace(0, 100, 3000), args=(mu,))
        ax.plot(sol.y[0], sol.y[1], linewidth=0.5, alpha=0.7)

    ax.set_xlabel('$x$')
    ax.set_ylabel('$\\dot{x}$')
    ax.set_title(f'Phase Portrait: $\\mu = {mu}$')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('van_der_pol.pdf', dpi=150, bbox_inches='tight')
plt.close()

# Calculate limit cycle amplitude
mu_test = 1.0
sol_lc = solve_ivp(van_der_pol, (0, 100), [0.1, 0], t_eval=np.linspace(0, 100, 5000), args=(mu_test,))
limit_cycle_amp = np.max(np.abs(sol_lc.y[0][-1000:]))
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{van_der_pol.pdf}
\caption{Van der Pol oscillator: time series (top) and phase portraits (bottom) showing limit cycles.}
\end{figure}

The limit cycle amplitude for $\mu = 1$ is approximately \py{f'{limit_cycle_amp:.3f}'}.

\chapter{Predator-Prey Dynamics}

\section{Lotka-Volterra Equations}
\begin{align}
\frac{dx}{dt} &= \alpha x - \beta xy \\
\frac{dy}{dt} &= \delta xy - \gamma y
\end{align}

\begin{pycode}
def lotka_volterra(t, state, alpha=1.0, beta=0.1, delta=0.075, gamma=1.5):
    x, y = state
    return [alpha*x - beta*x*y, delta*x*y - gamma*y]

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

t_span = (0, 50)
t_eval = np.linspace(0, 50, 1000)
sol = solve_ivp(lotka_volterra, t_span, [10, 5], t_eval=t_eval)

# Time series
ax = axes[0]
ax.plot(sol.t, sol.y[0], 'b-', label='Prey $x$')
ax.plot(sol.t, sol.y[1], 'r-', label='Predator $y$')
ax.set_xlabel('Time $t$')
ax.set_ylabel('Population')
ax.set_title('Population Dynamics')
ax.legend()
ax.grid(True, alpha=0.3)

# Phase portrait
ax = axes[1]
for x0, y0 in [(10, 5), (20, 10), (30, 5), (5, 2)]:
    sol_pp = solve_ivp(lotka_volterra, t_span, [x0, y0], t_eval=t_eval)
    ax.plot(sol_pp.y[0], sol_pp.y[1], linewidth=1)
ax.set_xlabel('Prey $x$')
ax.set_ylabel('Predator $y$')
ax.set_title('Phase Portrait')
ax.grid(True, alpha=0.3)

# Equilibrium
x_eq = 1.5 / 0.075
y_eq = 1.0 / 0.1
ax.plot(x_eq, y_eq, 'ko', markersize=8)

# Conservation quantity
ax = axes[2]
H = 0.075*sol.y[0] + 0.1*sol.y[1] - 1.5*np.log(sol.y[0]) - 1.0*np.log(sol.y[1])
ax.plot(sol.t, H, 'g-')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$H(x, y)$')
ax.set_title('Conserved Quantity')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lotka_volterra.pdf', dpi=150, bbox_inches='tight')
plt.close()

period_peaks = []
for i in range(1, len(sol.y[0])-1):
    if sol.y[0][i] > sol.y[0][i-1] and sol.y[0][i] > sol.y[0][i+1]:
        period_peaks.append(sol.t[i])
if len(period_peaks) > 1:
    lv_period = np.mean(np.diff(period_peaks))
else:
    lv_period = 0
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{lotka_volterra.pdf}
\caption{Lotka-Volterra predator-prey model: populations oscillate with period $T \approx \py{f'{lv_period:.2f}'}$.}
\end{figure}

\chapter{Numerical Results}

\begin{pycode}
# Summary of eigenvalue classifications
classifications = [
    ('Real, opposite signs', 'Saddle', 'Unstable'),
    ('Real, both negative', 'Stable node', 'Stable'),
    ('Real, both positive', 'Unstable node', 'Unstable'),
    ('Complex, $\\text{Re} < 0$', 'Stable spiral', 'Stable'),
    ('Complex, $\\text{Re} > 0$', 'Unstable spiral', 'Unstable'),
    ('Pure imaginary', 'Center', 'Neutral'),
]
\end{pycode}

\begin{table}[htbp]
\centering
\caption{Equilibrium classification by eigenvalues}
\begin{tabular}{@{}lll@{}}
\toprule
Eigenvalues & Type & Stability \\
\midrule
\py{classifications[0][0]} & \py{classifications[0][1]} & \py{classifications[0][2]} \\
\py{classifications[1][0]} & \py{classifications[1][1]} & \py{classifications[1][2]} \\
\py{classifications[2][0]} & \py{classifications[2][1]} & \py{classifications[2][2]} \\
\py{classifications[3][0]} & \py{classifications[3][1]} & \py{classifications[3][2]} \\
\py{classifications[4][0]} & \py{classifications[4][1]} & \py{classifications[4][2]} \\
\py{classifications[5][0]} & \py{classifications[5][1]} & \py{classifications[5][2]} \\
\bottomrule
\end{tabular}
\end{table}

\chapter{Conclusions}

\begin{enumerate}
    \item First-order ODEs: direction fields reveal solution behavior
    \item Second-order ODEs: damping ratio determines oscillation character
    \item Phase portraits: eigenvalues classify equilibrium types
    \item Limit cycles: nonlinear systems can have isolated periodic orbits
    \item Predator-prey: conservative systems show closed orbits
\end{enumerate}

\end{document}