Numerical ODEs: Euler, Runge-Kutta, Shooting & Systems

TAKEAWAYS:

• ODEs model real-world change: From population growth to bank interest, ODEs describe dynamic systems where rates depend on variables (e.g., $\frac{dy}{dx}=f(x,y)$).
• Euler’s method is intuitive but limited: It approximates solutions step-by-step with $y_{n+1}=y_n+h\cdot f(x_n,y_n)$, but errors accumulate for stiff equations.
• Runge-Kutta improves accuracy: The 4th-order method (RK4) balances precision and computation, using weighted averages of slopes.
• Shooting methods convert BVPs to IVPs: They guess initial conditions for boundary-value problems (BVPs) and refine iteratively.
• Systems of ODEs model interactions: Coupled equations (e.g., predator-prey models) require vectorized methods like Euler or RK4.
• Error control is critical: Local truncation error (LTE) and global error grow with step size $h$; adaptive methods adjust $h$ dynamically.

1. Ordinary Differential Equations (ODEs): Basics and Classification

ODEs describe how a quantity $y$ changes with respect to an independent variable $x$. They are classified by:

• Order: Highest derivative (e.g., $y^{\prime \prime }+y=0$ is 2nd-order).
• Linearity: Linear ODEs have $y$ and its derivatives to the first power (e.g., $y^{\prime }+p(x)y=q(x)$).
• Initial/Boundary Value Problems (IVP/BVP):
  • IVP: Specifies $y$ at a single point (e.g., $y(0)=1$).
  • BVP: Specifies $y$ at multiple points (e.g., $y(0)=1$, $y(1)=2$).

1st-order linear ODE: \( y' + p(x)y = q(x) \)
2nd-order linear ODE: \( y'' + a y' + b y = f(x) \)
Nonlinear ODE: \( y' = x^2 + y^2 \)
IVP example: \( y' = 3x^2 + 1 \), \( y(1) = 2 \)
BVP example: \( y'' + xy' - xy = 2x \), \( y(0)=1 \), \( y(2)=10 \)

Real-World Example 1: Bank Loan Interest (IVP)

A bank calculates loan repayment using the ODE for compound interest: $$\frac{dA}{dt}=rA$$ where $A$ is the loan amount, $r$ is the interest rate, and $t$ is time. For a loan of Rs. 100,000 at 5% annual interest ($r=0.05$), the exact solution is $A(t)=100000e^{0.05t}$. Numerically solving this with Euler’s method approximates the loan balance over time.

2. Euler’s Method: The Simplest Numerical Approach

Euler’s method approximates solutions by linearizing the ODE over small steps $h$: $$y_{n+1}=y_n+h\cdot f(x_n,y_n)$$ where $f(x,y)=\frac{dy}{dx}$.

Worked Example 1: Euler’s Method for $y^{\prime }=3x^2+1$, $y(1)=2$, $h=0.2$

Goal: Estimate $y(2)$. Steps:

1. Start at $x_0=1$, $y_0=2$.
2. Compute $f(x_0,y_0)=3(1{)}^2+1=4$.
3. Update: $$y_1=y_0+h\cdot f(x_0,y_0)=2+0.2\cdot 4=2.8$$
4. Repeat for $x_1=1.2$, $y_1=2.8$: $$f(x_1,y_1)=3(1.2{)}^2+1=5.12$$ $$y_2=2.8+0.2\cdot 5.12=3.824$$
5. Continue to $x_5=2$: $$y_5\approx 5.296$$

x: [1.0, 1.2, 1.4, 1.6, 1.8, 2.0]
y: [2.0, 2.8, 3.824, 5.072, 6.566, 5.296] (approximate)
Exact solution at x=2: y(2) ≈ 5.333 (from integrating \( y = x^3 + x + C \))
Error: |5.296 - 5.333| ≈ 0.037

Advantages/Disadvantages of Euler’s Method

3. Runge-Kutta Methods: Higher Accuracy

Euler’s method has local truncation error (LTE) $O(h^2)$. Runge-Kutta methods reduce this error by averaging slopes at intermediate points. The 4th-order Runge-Kutta (RK4) is widely used: $$\begin{array}{rl}{k}_1& =h\cdot f(x_n,y_n)\\ k_2& =h\cdot f\left(x_n+\frac{h}{2},y_n+\frac{k_1}{2}\right)\\ k_3& =h\cdot f\left(x_n+\frac{h}{2},y_n+\frac{k_2}{2}\right)\\ k_4& =h\cdot f(x_n+h,y_n+k_3)\\ y_{n+1}& =y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\end{array}$$

Worked Example 2: RK4 for $y^{\prime }=3x^2+1$, $y(1)=2$, $h=0.2$

Goal: Compare RK4 to Euler’s method. Steps:

1. At $x_0=1$, $y_0=2$: $$\begin{gathered} k_1=0.2\cdot (3(1{)}^2+1)=0.8 \\ {k}_2=0.2\cdot (3(1.1{)}^2+1)=0.882 \\ {k}_3=0.2\cdot (3(1.1{)}^2+1)=0.882(\text{Note: }{x}_n+h/2=1.1) \\ {k}_4=0.2\cdot (3(1.2{)}^2+1)=0.972 \\ {y}_1=2+\frac{1}{6}(0.8+2\cdot 0.882+2\cdot 0.882+0.972)=2.8827 \end{gathered}$$
2. Repeat for $x_1=1.2$, $y_1=2.8827$: $$y_2\approx 3.9366$$
3. Continue to $x_5=2$: $$y_5\approx 5.333(\text{Matches exact solution})$$

x: [1.0, 1.2, 1.4, 1.6, 1.8, 2.0]
Euler y: [2.0, 2.8, 3.824, 5.072, 6.566, 5.296]
RK4 y:   [2.0, 2.8827, 3.9366, 5.1606, 6.5556, 5.3330]
Exact y: [2.0, 2.88, 3.936, 5.16, 6.556, 5.333]

Comparison: Euler vs. RK4

4. Shooting Method for Boundary Value Problems (BVPs)

BVPs specify conditions at multiple points (e.g., $y(a)=\alpha$, $y(b)=\beta$). The shooting method converts a BVP into an IVP by guessing an initial slope and adjusting iteratively.

Worked Example 3: Shooting Method for $y^{\prime \prime }+x{y}^{\prime }-xy=2x$, $y(0)=1$, $y(2)=10$

Steps:

1. Rewrite as a system of 1st-order ODEs: $$y^{\prime }=v,v^{\prime }=-xy+xv+2x$$
2. Guess $v(0)=c$ (initial slope). Use Euler’s method to solve the system from $x=0$ to $x=2$.
3. Compare $y(2)$ to the target (10). Adjust $c$ using regula falsi or Newton’s method until $y(2)\approx 10$.

flowchart TD
    A["Guess \( v(0) = c_0 \)"] --> B["Solve IVP using Euler/RK4"]
    B --> C["Check \( y(2) \) vs. 10"]
    C -->|"If  y(2) < 10 "| D["Increase \( c \)"]
    C -->|"If  y(2) > 10 "| E["Decrease \( c \)"]
    D --> A
    E --> A

Real-World Example 2: Daraz Order Fulfillment (BVP)

Daraz’s warehouse uses BVPs to model inventory levels across regions. The ODE for inventory $I(t)$ at a warehouse satisfies: $$\frac{dI}{dt}=\text{Orders}-\text{Shipments}$$ with boundary conditions:

• $I(0)=\text{Initial stock}$,
• $I(T)=\text{Desired stock at time }T$. The shooting method helps optimize restocking schedules.

5. Systems of ODEs: Modeling Interacting Variables

Systems of ODEs describe coupled dynamics (e.g., predator-prey models, electrical circuits). The general form is: $$\frac{d\mathbf{y}}{dx}=\mathbf{f}(x,\mathbf{y}),\mathbf{y}(x_0)={\mathbf{y}}_0$$ where $\mathbf{y}=[y_1,y_2,\dots ,y_n{]}^T$.

Worked Example 4: Euler’s Method for a System

Solve: $$\frac{dy}{dx}=3x+\frac{y}{2},y(0)=1,h=0.1$$ Steps:

1. At $x_0=0$, $y_0=1$: $$\begin{gathered} f(x_0,y_0)=3(0)+\frac{1}{2}=0.5 \\ {y}_1=1+0.1\cdot 0.5=1.05 \end{gathered}$$
2. At $x_1=0.1$, $y_1=1.05$: $$\begin{gathered} f(x_1,y_1)=3(0.1)+\frac{1.05}{2}=0.525 \\ {y}_2=1.05+0.1\cdot 0.525=1.1025 \end{gathered}$$
3. Continue to $x=0.2$: $$y(0.2)\approx 1.1607$$

x: [0.0, 0.1, 0.2]
y: [1.0000, 1.0500, 1.1025, 1.1607]

Real-World Example 3: Kathmandu Traffic Flow (System of ODEs)

Traffic engineers model vehicle flow $Q(t)$ and density $\rho (t)$ using: $$\frac{d\rho }{dt}=-\frac{\rho Q}{\text{Length}}+\text{Entry rate}$$ $$\frac{dQ}{dt}=\text{Acceleration}-\text{Friction}$$ Euler’s method simulates congestion patterns to optimize signal timings.

6. Error Analysis and Step Size Control

• Local Truncation Error (LTE): Error per step (e.g., $O(h^2)$ for Euler).
• Global Error: Accumulated over all steps (e.g., $O(h)$ for Euler).
• Adaptive Methods: Adjust $h$ dynamically (e.g., Runge-Kutta-Fehlberg (RKF45)).

Step size (h): [0.1, 0.05, 0.01]
Global Error:  [0.05, 0.025, 0.005]

In the Real World

1. eSewa/Khalti Transactions (ODEs for Fraud Detection)

   • Banks use ODEs to model transaction rates $\frac{dT}{dt}=f(T,t)$, where $T$ is transaction volume. Anomalies (e.g., sudden spikes) trigger fraud alerts. Euler’s method approximates $T(t)$ to detect unusual patterns.

2. Pathao Driver Routing (Shooting Method for BVPs)

   • Pathao’s algorithm solves BVPs to find optimal routes between pickup/drop points. The ODE for distance $D(t)$ satisfies: $$\frac{d^{2}D}{d{t}^2}=\text{Acceleration}(t),D(0)=0,D(T)=\text{Destination}$$ The shooting method guesses initial acceleration to minimize travel time.

3. NTC Electricity Demand Forecasting (Systems of ODEs)

   • NTC models demand $D(t)$ and supply $S(t)$ using: $$\begin{gathered} \frac{dD}{dt}=\alpha D+\beta \cdot \text{Weather}(t) \\ \frac{dS}{dt}=\gamma S+\delta \cdot \text{Generation}(t) \end{gathered}$$ Euler’s method predicts blackout risks during peak hours.

Exam Tip

1. For IVPs:

   • Always show the Euler/RK4 update formula explicitly.
   • Compare numerical results to exact solutions (if possible) to highlight errors.
   • Example: In the exam, if asked to solve $y^{\prime }=x^2+y$ with Euler, write: $$y_{n+1}=y_n+h(x_n^2+y_n)$$ and compute at least 3 steps.

2. For BVPs:

   • Explain the shooting method’s iteration process (guess, solve IVP, adjust).
   • Use regula falsi or Newton’s method for adjustment (mention convergence criteria).

3. For Systems of ODEs:

   • Rewrite higher-order ODEs as 1st-order systems (e.g., $y^{\prime \prime }=f(x,y,y^{\prime })$ becomes $y^{\prime }=v$, $v^{\prime }=f(x,y,v)$).
   • Apply Euler/RK4 component-wise.

4. Common Pitfalls:

   • Forgetting to update $x$ in each step ($x_{n+1}=x_n+h$).
   • Misapplying boundary conditions in BVPs (e.g., mixing IVP/BVP setups).
   • Ignoring units in real-world examples (e.g., $h$ in years for population models).

5. High-Score Strategies:

   • Draw diagrams: Sketch the ODE’s slope field and numerical solution trajectory.
   • Label clearly: In RK4, write $k_1,k_2,\dots$ with their formulas.
   • Discuss limitations: For Euler, mention "large errors for stiff equations" or "requires small $h$".

Visual Summary of Methods:

mindmap
  root((Numerical ODEs))
    Euler
      Formula: y_{n+1} = y_n + h*f(x_n, y_n)
      Error: O(h)
      Use: Simple problems
    RK4
      Formula: Weighted avg of 4 slopes
      Error: O(h^4)
      Use: High accuracy needed
    Shooting Method
      Steps: Guess -> Solve IVP -> Adjust
      Use: BVPs
    Systems
      Form: dy/dx = f(x, y1, y2, ...)
      Method: Apply Euler/RK4 to each component