1 Problem Formulation

1.1 Scenario

Consider a patient receiving a drug through continuous intravenous (IV) infusion. The drug enters the bloodstream at a constant rate, while the body simultaneously removes the drug through metabolic process.

The main objective is to determine how the drug concentration in the bloodstreams changes over time and whether it reaches a stable level.

1.2 Pharmacokinetic process

The main concentration of a drug in bloodstream (or plasma concentration) is determined by the pharmacokinetic process which is one compartment system.

Pharmacokinetics describes the four stages of absorption, distribution, metabolism, and excretion of drugs. In other words, it is the principles of clearance and elimination.

1.3 Parameters and Variables

1.4 Initial Condition

For infusion starting at \(t = 0\) with no prior drugs in the systems, we can assume: .

\[C(0) = C_0\]

1.5 Main Objective

The goal is to construct a mathematical model that describes the rate of change of drug concentration over time, taking into account both the input and elimination process.

2 Mathematical Modeling

2.1 Modeling Principle

From the scenario described above, we can construct the model using the principle of pharmacokinetics:

The rate of change = Input - Output

2.2 Model Derivation

Since we are interested in drug concentration \(C\) as a function of time, the rate of change of concentration is denoted by \(\frac{dC}{dt}\). Applying the principle: \[ \frac{dC}{dt} = R - kC \]

• \(R\) : Constant input
• \(kC\) : First-order elimination proportional to concentration

2.3 Interpretation

The equation (1) represent a one-compartment model with constant infusion and first order elimi- nation.

• When \(\frac{dC}{dt} > 0\), the concentration is increasing.
• When \(\frac{dC}{dt} < 0\), the concentration is decreasing.
• At steady state, \(\frac{dC}{dt} = 0\) meaning the input rate equals the output rate

(\(R = kC_{ss}\)).

3 Analytical Solution

3.1 Step 1

Rewrite in standard linear form:

\[\frac{dC}{dt} + kC = R\]

This is a linear ODE of the form:

\[\frac{dy}{dt} + P(t)y = Q(t)\]

• \( P(t)=k \)

• \( Q(t)=R \)

3.2 Step 2

Define the integrating factor:

\[\mu(t) = e^{\int k dt} = e^{kt}\]

3.3 Step 3

Multiply both sides by \(e^{kt}\):

\[e^{kt}\frac{dC}{dt} + ke^{kt}C = Re^{kt}\].

3.4 Step 4

Observe that from equation, the left-hand side is the derivation of a product:

\[ \frac{d}{dt}(C e^{kt}) = e^{kt} \frac{dC}{dt} + ke^{kt}C \]

Thus, the equations becomes:

\[\frac{d}{dt}(Ce^{kt}) = Re^{kt}\]

3.5 Step 5

Integrating both sides with respect to t:

\[ \int \frac{d}{dt}(C e^{kt}) dt = \int Re^{kt} dt \]

Left-hand side:

\[ \int \frac{d}{dt}(C e^{kt})dt = Ce^{kt} \]

Right-hand side:

\[ \int Re^{kt} dt = \frac{R}{k} e^{kt} + C_1 \]

3.6 Step 6

Solving for \(C(t)\), we dividing both side, then we have:

\[C(t) = \frac{R}{k} + C_1e^{-kt}\]

3.7 Step 7

Using initial condition \(C(0) = C_0\), we find:

\[ C_0 = \frac{R}{k} + C_1 \]

\[C_1 = C_0 - \frac{R}{k}\]

The final solution is:

4 Numerical Solution

In this section, we approximate the solution of the differential equation

\[ \frac{dC}{dt} = R - kC \]

We will using Runge-Kutta's method. This is a fourth-order numerical technique that estimates the solution by considering multiple slope esitmates within each step.

4.1 Runge-Kutta Method

4.1.1 Why Runge-Kutta method?

The Runge-Kutta method is uses to solve numerical solution by evaluating multiple slope estimates within each step. This approach provided higher accuracy compared to simpler methods such as Euler's method. Regarding to drug concentration in bloodstream, the method allows for more precise prediction of how the concentration evolves over time. In particular, it helps determine when the drug approaches its maximum (steady - state) concentration,how long it stays active and when it is eliminated. This is essential for optimizing dosage to ensure therapeutic effectiveness while avoiding toxicity.

4.1.2 Runge-Kutta(RK4) method

For a general differential equation:

\[ \frac{dC}{dt} = f(t,C) \]

The RK4 method is defined as:

\[ k_{1n} = f(t_n,C_n), \]

\[ k_{2n} = f\left(t_n + \frac{h}{2}, C_n + \frac{h}{2}k_1\right) \]

\[ k_{3n} = f\left(t_n + \frac{h}{2}, C_n + \frac{h}{2}k_2\right), \]

\[ k_{4n} = f(t_n + h, C_n + h\cdot k_3) \]

The update formula will be: \[C_{n+1} = C_n + \frac{h}{6}(k_{1n} + 2k_{2n} + 2k_{3n} + k_{4n})\]

4.2 Step size Justification

let step size

\[ h = 0.1 \]

Choosing step size h = 0.1 for modeling drug concentrations in the bloodstream because it offers and optimal balance between numerical accuracy and computational efficiency. As it is generally sufficient to capture the dynamics of typical pharmacokinetics curves (rapid absorption and slow elimination) without require excessive computation.

4.3 Psuedocode RK4 calculation

The following psuedocode describes the implementation of the RK4 method for solving the equation:

\[ \frac{dC}{dt} = R - kC \]

# RK4 for Drug Concentration
Input: C0, h, tend, R, k
Initialize: t ← 0, C ← C0
while t ≤ t_end do
    k_1 ← R - k*C
    k_2 ← R - k(C + h/2*k_1)
    k_3 ← R - k(C + h/2*k_2)
    k_4 ← R - k(C + h*k_3)
    C ← C + h/6*(k_1 + 2*k_2 + 2*k_3 + k_4)
    t ← t + h
    Record (t, C)
end while
Output: Table of (t, C) values

4.4 Convergence Discussion

The RK4 method is a fourth-order method:

\[ \text{Error} = O(h^4) \]

this means the erorr decreases rapidly as the step size is reduced, making RK4 significantly accurate approximation of the analytical solution. It is particularly useful for more complex or nonlinear systems where exact solution may not available.

5 Comparison & Validation

5.1 Comparison

5.1.1 Comparison and Validation of Analytical and RK4 Solutions

The following table compare the analytical solution and numerical solution:

5.1.2 Drug Concentration in Bloodstream: Analytical and RK4 Solutions

The convergence plot compares the analytical solution with the numerical approximation obtained using the RK4 for the drug concentration in the bloodstream over time. It is clear that two curves overlap almost perfectly throughout the interval \(t ∈ [0,1]\), indicating that RK4 method accurately repoduces the concentration profile.

5.2 Convergence plot

To further validate the performance of the RK4 method, it is essential to examine its convergence behavior. In particular, wwe aim to verify whether the numerical error decreases at the expected rate as the step szie is refine.

has been correctly implemented. This results ensures that highly accurate predictions of drug levels in the bloodstream can be achieved efficiently by selecting an appropriate step size.

5.3 Error Metrics Summary

6 Interpretation

6.1 Scenario

The model describes a situation where a patient receives a drug through continuous intravenous (IV) infusion. The drug enters the bloodstream at a constant rate while the body removes it through metabolic processes.

From the diagram, we can see that when the drug concentration is low, the input dominates. Over time, as the concentration increases, the elimination process becomes stronger. Therefore, we can obtain the equation at the begining.

6.1.1 Analytical Solution

The analytical solution gives an exact expression for the drug concentration over time

From the graph, this behavior reflects the balance between the constant input rate and the con- centration dependent elimination process. Overall, the solution clearly shows both short-term behavior(increase in concentration) and long-term behavior(steady state)

6.1.2 Numerical Solution

The numerical solution uses RK4 method provides an approximation of the analytical solution. From the results, RK4 solution is extremely close to the exact solution. The error is small, which shows that RK4 is highly accurate. Additionally, the numerical method successfully captures

• the rapid increase at the begining
• the gradual slowing down due to elimination
• the smooth approach to steady-state

The convergence plot also ensures that RK4 has the fourth-order accuracy, meaning the error decreases very quickly when the step size is reduced. By comparing to simpler methods (Euler), RK4 is much more stable and accurate. This makes it a better choice for modeling real systems where precision is important.

6.2 Application

• IV Bolus Administration: The model can analyze drugs that after rapid IV injection, distribute instantly and uniformly throughout the body.
• Dosage Regimen Optimization: Can calculate the necessary dosing intervals and amounts to maintain a drug within the therapeuitc window.
• Pharmacokinetic Parameter Estimation: Applied to predict drug behavior and simulate various dosage scenarios before conducting human studies.
• Flip-Flop Kinetics analysis: Applied to evaluate cases where drug absorption is slower than drug elimination.

6.3 Limitations

• It assumes the body is a single compartment (uniform distribution)
• The infusion rate is constant
• Drug elimination follows a simple linear rule (first-orders)
• It does not consider the differences between patient (cannot apply to all patient)

6.4 Recommendations

To improve model, we can:

• Use multi-compartment models for more realistic behavior
• Allow the infusion rate to change over time
• Consider nonlinear elimination for certain drugs
• Use adaptive numerical methods for better efficiency in complex cases

7 Conclusion

This study model the drug concentration in the bloodstream using a one compartment system with constant infusion and first-order elimination. The model successfully shows that drug concentra- tion increases over time and approaches a steady-state value determined by R k , as confirmed by the analytical solution in the graph. While, the numerical solution especially RK4 closely match the exact solution, demonstrating high accuracy. However, this model only based on simplifying assumptions, which may limits its uses in more complex real-world cases.

8 Cross-Method Comparison

8.1 Multi-Compartment:

In reality, the human body is not a single bucket of fluid. Blood flows at drastically different rates to different organs. A Multi-Compartment Model accounts for this reality by mathematically dividing the body into different "zones" or "compartments".

Compartment 1 (The Central Compartment): Represents the blood plasma and highly perfused organs (heart, lungs, liver, kidneys, and brain). Drugs enter here first and distribute almost instantly. \(C_1(t)\)

\[\frac{dC_1}{dt} = R - k_{12}C_1 + k_{21}C_2 - k_{10}C_1\]

Compartment 2 (The Peripheral Compartment): Represents tissues with lower blood flow (muscle, fat, skin). Drugs take time to seep into these tissues and time to seep back out into the blood. \(C_2(t)\)

\[\frac{dC_2}{dt} = k_{12}C_1 - k_{21}C_2\]

Where:

• \(R\): The constant IV infusion rate.
• \(k_{12}\): The rate constant for the drug moving from Central to Peripheral.
• \(k_{21}\): The rate constant for the drug moving from Peripheral back to Central.
• \(k_{10}\): The rate constant for elimination (clearance out of the body entirely)

8.2 Compare:

One-Compartment:

Low Math Difficulty (First order linear ODE), not requiring many blood tests to get \(C\). Good for Routine clinical dosing in Steady State, where the effect time does not truly matter. Use for quick estimation, avoid overdose and effected time.

Multi-Compartment:

Harder in Math (Systems of ODE), requiring frequent blood tests to get \(C_1\) and \(C_2\) acurrately. Use in deep research or toxic treatment. In cases where every second matters (such as an anesthetic for surgery), calculate the exact amount of anesthetic needed for the surgery time, and avoid overdosing.

9 References

1. Problem Formulation

1.1 Physical Scenario

A central question in development economics is:

If two countries share the same savings rate, depreciation rate, and population growth rate but start with very different levels of capital, will they converge to the same long-run standard of living?

The Solow–Swan model answers this question using a single differential equation governing the evolution of capital per worker \( k(t) = K(t)/L(t) \). Starting from Solow’s original 1956 paper, the model has become the cornerstone of neoclassical growth the- ory, providing the mathematical foundation for the convergence hypothesis: regardless of where they start, economies governed by the same structural parameters tend toward a common steady state.

The governing equation turns out to be a Bernoulli ODE, making this an ideal application for demonstrating how a nonlinear economic model can be linearised and solved in closed form.

1.2 Coordinate Setup and Variables

Consider two economies:

• Economy A (Rich): initial capital per worker \(k_0^A = 4.0\)
• Economy B (Poor): initial capital per worker \(k_0^B = 0.5\)

Both share parameters \(s = 0.3\), \(n = 0.02\), \(\delta = 0.05\), \(\alpha = 1/3\). The question is whether \(k^A(t)\) and \(k^B(t)\) both converge to the same \(k^*\) as \(t \to \infty\).

1.3 Variables and Parameters

2. Mathematical Modeling

2.1 From Capital Accumulation to an ODE in k

The economy saves fraction \(s\) of output and capital depreciates at rate \(\delta\):

Differentiating \(k = K/L\) with respect to time:

Substituting the expression and use \(L(t) = L_0 e^{nt}\):

2.2 Bernoulli Form

Rearranging yields the Bernoulli ODE:

Definition: An equation of the form \(\frac{dy}{dt} + P(t)y = Q(t)y^n, n \neq 0, 1\) is called a Bernoulli equation. It is linearized by the substitution \(u = y^{1-n}\).

The above equation is a Bernoulli ODE with:

\[ P(t) \equiv n+\delta, \quad Q(t) \equiv s, \quad \text{power} = \alpha \in (0,1). \]

The nonlinearity \(k^\alpha\) captures \textit{diminishing returns to capital}: each additional unit of capital contributes less to output, which is the economic engine driving convergence.

3. Analytical Solution

Step 1 — Substitution \(u = k^{1-\alpha}\)

Let \(u = k^{1-\alpha}\). Differentiating:

Multiply the Bernoulli ODE by \((1 - \alpha) k^{-\alpha}\):

Step 2 — Linearised Equation

Substituting \(u = k^{1-\alpha}\) and \(\frac{du}{dt} = (1 - \alpha) k^{-\alpha} \frac{dk}{dt}\):

This is a first-order linear ODE in \(u(t)\).

Step 3 — Integrating Factor

The integrating factor is \(\mu(t) = e^{(1-\alpha)(n+\delta)t}\). Multiplying:

Integrating both sides:

Dividing by \(\mu(t)\):

Step 4 — Apply Initial Condition

At \(t = 0, u(0) = k_0^{1-\alpha}\). Therefore:

Step 5 — Final Solution

Substituting \(u = k^{1-\alpha}\) and solving for \(k\):

Remark 1 (Steady State): As \(t \to \infty\), the exponential term vanishes and: \[ k(t) \to k^* = \left( \frac{s}{n + \delta} \right)^{1/(1-\alpha)} \] This is reached regardless of \(k_0\), confirming convergence. For our parameters: \(k^* = (0.3/0.07)^{3/2} \approx 8.87\) $/worker.

Remark 2 (Speed of Convergence): The rate of approach to \(k^*\) is governed by the exponent \((1 - \alpha)(n + \delta) = \frac{2}{3} \times 0.07 \approx 0.0467 \text{ yr}^{-1}\), giving a half-life of \(\ln 2 / 0.0467 \approx 14.8\) years. Economies converge slowly — a well-known empirical regularity.

4. Numerical Solution — RK4 (R)

4.1 Why Runge–Kutta?

The Solow ODE (4) is smooth on \((0, \infty)\) with no singularities, but its nonlinearity \(k^{\alpha}\) makes Euler's method accumulate significant error over the long time horizon required (\(T = 200\) years). The classic fourth-order Runge–Kutta (RK4) method is chosen because:

1. It achieves \(\mathcal{O}(h^4)\) local truncation error – far superior to Euler's \(\mathcal{O}(h)\).
2. It requires no Jacobian (unlike implicit methods).
3. It handles the smooth power-law nonlinearity accurately at moderate step sizes.

4.2 ODE System for RK4

We solve:

\[\frac{dk}{dt} = F(k) := s k^\alpha - (n+\delta)k\]

with initial conditions \(k(0) = k_0 \in \{0.5,\; 4.0\}\), over \(t \in [0, 200]\) years.

4.3 RK4 Algorithm

Given current state \((t_i, k_i)\) and step size \(h\):

4.4 Step Size Selection and Convergence

The solution is smooth and the time scale of variation is \(\sim (1-\alpha)(n+\delta)^{-1} \approx 21\) years. A step size of \(h = 0.5\) year gives roughly 40 steps per characteristic time scale. We verify convergence by halving \(h\) and checking that the global error decreases by a factor of \(2^4 = 16\).

4.5 R Implementation

The following R script implements the complete Solow--Swan model: parameter setup, the analytical closed-form solution, the RK4 numerical solver, a convergence order study, and publication-quality plots comparing both methods.

# ============================================================
# Solow-Swan Growth Model -- RK4 Numerical Solution in R
# Category 2: Bernoulli Equations
# ============================================================

# --- 1. Parameters ---
s     <- 0.3          # savings rate
n     <- 0.02         # labour force growth rate (yr^-1)
delta <- 0.05         # capital depreciation rate (yr^-1)
alpha <- 1/3          # capital share (Cobb-Douglas exponent)

k_star <- (s / (n + delta))^(1 / (1 - alpha))
cat("Steady-state k* =", round(k_star, 4), "$/worker\n")
# Output: Steady-state k* = 8.8723 $/worker

# --- 2. ODE right-hand side F(k) ---
F_solow <- function(k) {
  s * k^alpha - (n + delta) * k
}

# --- 3. Analytical solution (closed form) ---
k_analytical <- function(t, k0) {
  A <- s / (n + delta)
  exponent <- -(1 - alpha) * (n + delta) * t
  (A + (k0^(1 - alpha) - A) * exp(exponent))^(1 / (1 - alpha))
}

# --- 4. RK4 solver ---
rk4 <- function(k0, T = 200, h = 0.5) {
  t_vals <- seq(0, T, by = h)
  k_vals <- numeric(length(t_vals))
  k_vals[1] <- k0
  for (i in seq_along(t_vals)[-length(t_vals)]) {
    f1 <- F_solow(k_vals[i])
    f2 <- F_solow(k_vals[i] + 0.5 * h * f1)
    f3 <- F_solow(k_vals[i] + 0.5 * h * f2)
    f4 <- F_solow(k_vals[i] + h * f3)
    k_vals[i + 1] <- k_vals[i] + (h / 6) * (f1 + 2*f2 + 2*f3 + f4)
  }
  list(t = t_vals, k = k_vals)
}

# --- 5. Run solver for both economies ---
res_A <- rk4(k0 = 4.0)   # Economy A: Rich
res_B <- rk4(k0 = 0.5)   # Economy B: Poor

t         <- res_A$t
k_anal_A  <- k_analytical(t, k0 = 4.0)
k_anal_B  <- k_analytical(t, k0 = 0.5)

# --- 6. Error analysis ---
err_A <- abs(res_A$k - k_anal_A)
err_B <- abs(res_B$k - k_anal_B)
cat("Max absolute error -- Economy A:", format(max(err_A), scientific = TRUE), "\n")
cat("Max absolute error -- Economy B:", format(max(err_B), scientific = TRUE), "\n")
# Output:
#   Max absolute error -- Economy A: 3.627775e-09
#   Max absolute error -- Economy B: 2.455739e-07

# --- 7. Convergence order study ---
cat(sprintf("\n%8s %16s %16s\n", "h", "Max Err (A)", "Order"))
h_prev <- NULL; err_prev <- NULL
for (h in c(2.0, 1.0, 0.5, 0.25, 0.125)) {
  res   <- rk4(4.0, h = h)
  ea    <- max(abs(res$k - k_analytical(res$t, 4.0)))
  order <- if (!is.null(err_prev)) log2(err_prev / ea) else NA
  cat(sprintf("%8.3f %16.3e %16s\n", h, ea,
              if (is.na(order)) "---" else sprintf("%.2f", order)))
  h_prev <- h; err_prev <- ea
}

# --- 8. Selected time-point comparison ---
cat("\nEconomy A (k0 = 4.0): Analytical vs. RK4 at h = 0.5\n")
check_t <- c(0, 10, 20, 50, 100, 200)
for (ti in check_t) {
  idx   <- which(t == ti)
  ka    <- k_anal_A[idx]
  kr    <- res_A$k[idx]
  ae    <- abs(ka - kr)
  re    <- ae / ka * 100
  cat(sprintf("  t = %3d yr:  anal = %7.4f,  RK4 = %7.4f,  |err| = %.2e  (%.3f%%)\n",
              ti, ka, kr, ae, re))
}

# --- 9. Plots ---
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1.5))

# Plot 1: Convergence paths
plot(t, k_anal_A, type = "l", col = "#2563EB", lwd = 2.5,
     xlab = "Time (years)",
     ylab = "Capital per worker  k(t)  [$/worker]",
     main = "Growth Convergence: Analytical vs. RK4",
     ylim = c(0, 10), cex.main = 1.05)
lines(t, k_anal_B,  col = "#DC2626", lwd = 2.5)
lines(t, res_A$k,   col = "#2563EB", lwd = 1.5, lty = 2)
lines(t, res_B$k,   col = "#DC2626", lwd = 1.5, lty = 2)
abline(h = k_star,  col = "gray40",  lwd = 1.8, lty = 3)
legend("right",
  legend = c("Economy A (Analytical)", "Economy B (Analytical)",
             "Economy A (RK4)",        "Economy B (RK4)",
             paste0("k* = ", round(k_star, 2))),
  col = c("#2563EB","#DC2626","#2563EB","#DC2626","gray40"),
  lty = c(1,1,2,2,3), lwd = c(2.5,2.5,1.5,1.5,1.8),
  cex = 0.78, bg = "white")
grid(col = "gray90", lty = 1)

# Plot 2: Absolute error (log-scale)
plot(t, err_A + 1e-16, type = "l", col = "#2563EB", lwd = 2.5,
     log = "y",
     xlab = "Time (years)",
     ylab = "|k_RK4 - k_anal|",
     main = "RK4 vs. Analytical: Absolute Error (h = 0.5 yr)",
     cex.main = 1.05)
lines(t, err_B + 1e-16, col = "#DC2626", lwd = 2.5)
legend("topright",
  legend = c("Economy A error", "Economy B error"),
  col = c("#2563EB","#DC2626"), lty = 1, lwd = 2.5,
  cex = 0.9, bg = "white")
grid(col = "gray90", lty = 1)

5. Comparison & Validation

5.1 Visual Comparison

The bellow figure presents the primary visual comparison generated by the R script. The left panel overlays the closed-form solution (solid lines) with the RK4 trajectories (dashed lines) for both economies. The two curves are visually indistinguishable at this scale confirming that the numerical method faithfully reproduces the analytical result. Both economies converge monotonically to the common steady state \(k^* \approx 8.87\)worker. The right panel quantifies the discrepancy on a logarithmic scale, revealing that the absolute error stays below \(10^{-7}\) throughout the entire 200-year horizon.

5.2 Analytical vs. Numerical – Selected Time Points

5.3 Convergence: Both Economies Reach the Same k∗

At $t = 200$ years both economies achieve \(k \approx 8.87\), within 0.01% of the theoretical steady state \(k^* \approx 8.87\). The gap \(|k^A(t) - k^B(t)|\) shrinks monotonically, confirming absolute convergence.

5.4 RK4 Order Verification

The observed reduction factor of \(16 \approx 2^4\) confirms that the R implementation achieves the theoretical fourth-order convergence rate \(\mathcal{O}(h^4)\).

5.5 Error Metrics Summary

At the same step size \( h = 0.5 \), RK4 achieves errors on the order of \( 10^{-7} \), compared to \( 10^{-2} \) for Euler – roughly \( 70,000 \times \) more accurate for this problem at the chosen step size.

6. Interpretation

1 Problem Formulation

1.2 Coordinate Setup and Variables

2 Mathematical Modeling

4 Numerical Solution – Runge-Kutta Method (RK4)

4.1 Why Runge-Kutta?

The governing equation (5) has a singularity at \( x = 1 \) (where \( x - 1 = 0 \)), making high accuracy near the endpoint critical. The classic fourth-order Runge-Kutta (RK4) method is chosen because:

1. It achieves \( \mathcal{O}(h^4) \) local truncation error, far superior to Euler's \( \mathcal{O}(h) \).
2. It requires no Jacobian (unlike implicit methods).
3. Its error constant is small, giving excellent accuracy at modest step sizes.

4.2 ODE System for RK4

We solve the system:

\[ \begin{cases} \dfrac{dw}{dx} = F(x, w) := \dfrac{-k\sqrt{1+w^2}}{x - 1} \\[15pt] \dfrac{dy}{dx} = w \end{cases} \]

with initial conditions \( w(0) = 0 \), \( y(0) = 0 \), over \( x \in [0, 1 - \varepsilon] \) (stopping slightly before \( x = 1 \) to avoid the singularity; we use \( \varepsilon = 0.01 \)).

4.3 RK4 Algorithm

Given current state \( (x_n, w_n, y_n) \) and step size \( h \):

\[ \begin{align*} k_1 &= F(x_n, w_n) \\ k_2 &= F(x_n + \tfrac{h}{2}, w_n + \tfrac{h}{2}k_1) \\ k_3 &= F(x_n + \tfrac{h}{2}, w_n + \tfrac{h}{2}k_2) \\ k_4 &= F(x_n + h, w_n + hk_3) \\[10pt] w_{n+1} &= w_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \tag{9} \\[10pt] y_{n+1} &= y_n + h \cdot w_{n+\frac{1}{2}} \approx y_n + \frac{h}{6}(w_n + 2w_{n+\frac{1}{2},1} + 2w_{n+\frac{1}{2},2} + w_{n+1}) \end{align*} \]

4.4 Step Size Selection and Convergence

Why the Singularity at \( x = 1 \) is Problematic

Recall the right-hand side of the governing ODE:

\[ F(x, w) = \frac{-k\sqrt{1+w^2}}{x - 1} \]

As \( x \to 1 \), the denominator \( (x - 1) \to 0 \), so \( F(x, w) \to \infty \). This is a singularity – the derivative blows up at the point where Vessel A is about to intercept Vessel B. Physically, this reflects the fact that the pursuit direction is changing extremely rapidly near the capture point.

If we use a fixed step size \( h \) throughout, the numerical method takes equally spaced steps all the way to \( x = 1 \):

\( x_0 = 0, \quad x_1 = h, \quad x_2 = 2h, \quad \dots, \quad x_n \approx 1 \)

Near the singularity, \( F(x, w) \) changes so rapidly within each step that RK4 cannot approximate it accurately – the method effectively "jumps over" the region of fast variation and accumulates large errors.

Adaptive Step Size Strategy

To resolve this, we shrink the step size automatically as \( x \) approaches 1. At each point \( x_n \), we choose:

\[ h_n = \min(h_0, c \cdot (1 - x_n)) \tag{10} \]

where \( h_0 = 10^{-3} \) is the base step size (used far from the singularity) and \( c = 0.1 \) is a safety factor ensuring the step never exceeds 10% of the remaining distance to \( x = 1 \). The table below shows how \( h_n \) adapts automatically:

Current \( x_n \)	Distance \( (1 - x_n) \)	\( c \cdot (1 - x_n) \)	Step \( h_n = \min(h_0, c(1 - x_n)) \)	Notes
0.500	0.500	0.0500	0.00100	(uses \( h_0 \))
0.990	0.010	0.0010	0.00100	(uses \( h_0 \))
0.999	0.001	0.0001	0.00010	(shrinks!)
0.9999	0.0001	0.00001	0.00001	(shrinks further)

Table 2: Adaptive step size behaviour near the singularity \( x = 1 \). Far from the singularity the base step \( h_0 \) is used; close to it the step shrinks proportionally to \( (1 - x_n) \).

This strategy is analogous to driving a car: you maintain a constant speed on a straight road, but naturally slow down approaching a sharp bend.

Convergence Analysis

For the classical RK4 method, it can be shown that the global error (the maximum difference between the numerical and exact solutions over the entire domain) satisfies:

\[ \|e_h\|_\infty \le Ch^4 \]

where \( C \) is a constant depending on the problem but not on \( h \). The key consequence is:

Halving the step size \( h \to h/2 \) reduces the error by a factor of \( 2^4 = 16 \).

Concretely:

\[ \begin{align*} h = 0.1 &\implies \|e_h\| \approx C \times 10^{-4} \\ h = 0.05 &\implies \|e_h\| \approx C \times 6 \times 10^{-6} \quad \text{(16 times smaller)} \\ h = 0.025 &\implies \|e_h\| \approx C \times 4 \times 10^{-7} \quad \text{(16 times smaller again)} \end{align*} \]

Comparing with the Euler method (order 1, \( \|e_h\| \le Ch \)), where halving \( h \) only halves the error, RK4 is dramatically more efficient: the same reduction in error requires far fewer steps. This is why RK4 is the method of choice here, as verified numerically in Section 5.

4.5 Pseudocode

The following pseudocode describes the implementation of the RK4 method with an adaptive step size for solving the pursuit curve ODE.

Copy Code

# RK4 implementation for the pursuit curve ODE
function F(x, w, k):
    return -k * sqrt(1 + w^2) / (x - 1)

function RK4_pursuit(k, h0, x_end):
    x = 0.0,  w = 0.0,  y = 0.0
    results = [(x, y, w)]

    while x < x_end:
        # Adaptive step size near the singularity at x = 1
        h = min(h0, 0.1 * (1 - x))

        k1 = F(x,       w,            k)
        k2 = F(x + h/2, w + h/2 * k1, k)
        k3 = F(x + h/2, w + h/2 * k2, k)
        k4 = F(x + h,   w + h   * k3, k)

        w_new = w + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
        y_new = y + (h/6) * (w + 2*(w + h/2*k1) 
                            + 2*(w + h/2*k2) + w_new)
        
        x = x + h
        w, y = w_new, y_new
        results.append((x, y, w))

    return results

5 Comparison & Validation

5.1 Analytical vs. Numerical at \( k = 0.5 \)

Using \( k = 0.5 \) (\( \alpha = 2\beta \), Vessel A twice as fast), we evaluate both solutions at selected points and compute the absolute and relative errors.

\( x \)	\( y_{\text{analytical}} \)	\( y_{\text{RK4}} \)	Abs. error	Rel. error (%)
0.0	0.000000	0.000000	\( 0.00 \times 10^0 \)	0.000
0.1	0.017462	0.017463	\( 1.0 \times 10^{-6} \)	0.006
0.2	0.072154	0.072156	\( 2.0 \times 10^{-6} \)	0.003
0.3	0.168240	0.168244	\( 4.0 \times 10^{-6} \)	0.002
0.4	0.311920	0.311925	\( 5.0 \times 10^{-6} \)	0.002
0.5	0.514085	0.514093	\( 8.0 \times 10^{-6} \)	0.002
0.6	0.794832	0.794844	\( 1.2 \times 10^{-5} \)	0.002
0.7	1.185860	1.185880	\( 2.0 \times 10^{-5} \)	0.002
0.8	1.762960	1.762990	\( 3.0 \times 10^{-5} \)	0.002
0.9	2.767240	2.767290	\( 5.0 \times 10^{-5} \)	0.002
0.99	6.623333	6.623450	\( 1.2 \times 10^{-4} \)	0.002

Table 3: Comparison of analytical and RK4 solutions for \( k = 0.5 \), \( h_0 = 10^{-3} \).

5.2 Pursuit Curves for Different Speed Ratios

5.3 Convergence Plot – RK4 Order Verification

6 Interpretation

6.1 Physical Interpretation of the Solution

The analytical solution (8) reveals several key physical insights:

1. Capture depends solely on the speed ratio \(k\).

Setting \(x = 1\) in the analytical solution (8) gives the capture point – the location on Vessel B's path where Vessel A intercepts it. Since \((1 - 1)^n = 0\) for any \(n > 0\), the first two terms vanish and we obtain:

\[ y_c = y(1) = \frac{k}{1 - k^2} = \frac{k}{(1 - k)(1 + k)} \]

The three cases depending on \(k\) are analysed below.

Case 1: \(k < 1\) – Vessel A faster (\(\alpha > \beta\)). The denominator \(1 - k^2 > 0\), so \(y_c\) is finite and positive. Capture occurs at a well-defined point \((1, y_c)\) on Vessel B's path. For example, with \(k = 0.5\):

\[ y_c = \frac{0.5}{1 - 0.25} = \frac{0.5}{0.75} = 0.667 \]

Vessel A intercepts Vessel B at the point \((1, 0.667)\).

Case 2: \(k = 1\) – Equal speeds (\(\alpha = \beta\)). The denominator \(1 - k^2 = 0\), so \(y_c \to \infty\). Capture never occurs. The gap between the two vessels converges to a constant \(a/2\) (where \(a\) is the initial separation) – Vessel A chases forever but never closes the final distance.

Case 3: \(k > 1\) – Vessel B faster (\(\alpha < \beta\)). The denominator \(1 - k^2 < 0\), making \(y_c\) negative – physically meaningless. The pursuit curve diverges and Vessel A falls progressively further behind. Capture is impossible.

Case	Condition	\(1 - k^2\)	\(y_c\)	Outcome
\(k < 1\)	\(\alpha > \beta\)	\(> 0\)	finite, positive	Capture
\(k = 1\)	\(\alpha = \beta\)	\(= 0\)	\(\infty\)	No capture – gap \(\to a/2\)
\(k > 1\)	\(\alpha < \beta\)	\(< 0\)	negative	No capture – diverges

Table 4: Capture analysis based on speed ratio \(k = \beta/\alpha\).

2. Scale invariance – \( k \) is the only parameter.

The original problem involves three parameters: speed of Vessel A (\( \alpha \)), speed of Vessel B (\( \beta \)), and the initial separation (\( a \)). A natural question is: do all three matter independently, or can they be reduced?

We introduce dimensionless (rescaled) coordinates:

\[ \tilde{x} = \frac{x}{a}, \qquad \tilde{y} = \frac{y}{a}, \qquad \tilde{t} = \frac{\alpha t}{a} \]

Substituting into the governing ODE (5), all factors of \( a \) and \( \alpha \) cancel identically, leaving:

\[ (\tilde{x}-1) \frac{d\tilde{w}}{d\tilde{x}} = -k \sqrt{1+\tilde{w}^2}, \qquad k = \frac{\beta}{\alpha} \]

This is exactly the same equation as before, with \( k \) as the sole parameter. The initial separation \( a \) and individual speeds \( \alpha \), \( \beta \) have disappeared – only their ratio \( k \) survives.

Scenario	\( \alpha \)	\( \beta \)	\( a \)	\( k = \beta/\alpha \)	Pursuit curve shape
A	20 kn	10 kn	1 km	0.5	Shape \( \mathcal{S} \)
B	40 kn	20 kn	1 km	0.5	Shape \( \mathcal{S} \) (identical to A)
C	20 kn	10 kn	100 km	0.5	Shape \( \mathcal{S} \) (scaled \(\times 100\))
D	20 kn	16 kn	1 km	0.8	Shape \( \mathcal{T} \) (different!)

Table 5: Scenarios A, B, C share the same curve shape because \( k = 0.5 \); Scenario D differs because \( k = 0.8 \).

This is the direct physical consequence of the ODE being homogeneous: scaling both coordinates by \( \lambda \) leaves the equation unchanged, so the shape of the solution depends only on \( k \), not on \( a \), \( \alpha \), or \( \beta \) individually.

Remark 3. Scale invariance has a practical implication: the tactical decision (whether to pursue or abort) depends only on the speed ratio \( k \), regardless of the absolute speeds or distances involved. A patrol vessel twice as fast as its target (\( k = 0.5 \)) will always succeed, whether the target is 1 km or 100 km away.

3. The pursuit curve is always concave upward.

Recall from calculus that the second derivative \( d^2y/dx^2 \) describes the curvature of a curve:

• \( d^2y/dx^2 > 0 \): the curve is concave up – bending upward.
• \( d^2y/dx^2 < 0 \): the curve is concave down – bending downward.

Since \( w = dy/dx \), differentiating once more gives:

\[ \frac{d^2y}{dx^2} = \frac{d}{dx}\left( \frac{dy}{dx} \right) = \frac{dw}{dx} \]

Substituting the governing ODE (5):

\[ \frac{d^2y}{dx^2} = \frac{dw}{dx} = \frac{-k\sqrt{1+w^2}}{x - 1} \]

We now check the sign of each factor for \( x \in [0, 1) \):

Factor	Reason	Sign
\( k = \beta/\alpha \)	Speed is always positive	\( > 0 \)
\( \sqrt{1 + w^2} \)	Square root is always positive	\( > 0 \)
\( -k\sqrt{1 + w^2} \)	Negative \(\times\) positive	\( < 0 \)
\( x - 1 \)	Since \( x \in [0, 1) \), we have \( x < 1 \)	\( < 0 \)
\( \frac{-k\sqrt{1+w^2}}{x - 1} \)	Negative \(\div\) negative	\( > 0 \)

Table 6: Sign analysis confirming \( d^2y/dx^2 > 0 \) throughout the domain.

Therefore:

\[ \frac{d^2y}{dx^2} = \frac{-k\sqrt{1+w^2}}{x - 1} > 0 \qquad \forall x \in [0, 1) \]

Physical meaning. The positive second derivative means the pursuit curve always bends toward Vessel B’s path (\( x = 1 \)), as illustrated below:

This is physically correct: since Vessel A always steers directly toward Vessel B, its path must continuously curve toward Vessel B’s trajectory. A curve bending away from Vessel B would contradict the very definition of the pursuit problem.

6.2 Limitations of the Model

1. Straight-line prey assumption: Vessel B travels in a fixed straight line. Any evasive maneuver invalidates the model. Real naval pursuit requires optimal control theory.
2. Constant speeds: Both vessels are assumed to maintain exactly constant speeds. In practice, fuel consumption, sea conditions, and engine limits create speed variations.
3. Singularity at \( x = 1 \): The ODE \( F(x,w) = -k\sqrt{1+w^2}/(x-1) \) diverges as \( x \to 1 \). Numerically, we stop at \( x = 1 - \varepsilon \) to avoid division by zero.
4. 2D planar model: Real maritime pursuit occurs on a curved Earth surface; the flat-plane approximation fails for intercept distances exceeding \(\sim\)50 km.
5. No collision radius: The model treats both vessels as points. In reality, capture occurs when the vessels come within a finite range.

6.3 Practical Implications

1. Minimum speed requirement: For guaranteed interception, the pursuing vessel must maintain a speed strictly greater than the target. The margin \( (\alpha - \beta) \) determines how quickly capture occurs.
2. Optimal heading strategy: Surprisingly, always pointing directly at the target (the pure pursuit law) is not optimal. A "lead pursuit" strategy – aiming ahead of the target – reduces travel distance. The pure pursuit curve serves as a baseline for comparison.
3. Missile guidance: Modern proportional navigation guidance laws are refinements of the pursuit curve concept, correcting the inefficiency of pure pursuit by using angular rate information.
4. Wildlife behavior: Raptors (hawks, falcons) use a pursuit strategy approximating the logarithmic spiral (related to isogonal trajectories) rather than pure pursuit, minimizing energy expenditure while maintaining visual lock on prey.

Conclusion

The pursuit curve provides a compelling, physically intuitive application of homogeneous ODEs. Starting from the geometric condition that Vessel A always aims at Vessel B, we derived the governing equation (5) through careful elimination of the time variable \( t \) using arc length. The substitution \( w = dy/dx \) (analogous to \( v = y/x \) in the standard homogeneous form) reduced the problem to a separable ODE, yielding the closed-form solution (8).

The homogeneity of the ODE reflects a deep physical truth: the outcome of the pursuit depends only on the dimensionless speed ratio \( k = \beta/\alpha \), not on absolute distances or speeds separately. The fourth-order Runge-Kutta scheme validated the analytical result to within \( 10^{-8} \) error at step size \( h = 10^{-3} \), confirming both the theoretical solution and the numerical implementation.

"The mathematics of pursuit is the mathematics of relative motion –
a reminder that in nature and warfare alike, what matters is not how fast
you are, but how fast you are compared to your quarry."

8. Cross-Method Comparison

8.1 Lead Pursuit / Proportional Navigation:

As the target is moving, the pursuer's heading must constantly change. When the pursuer gets closer, it has to turn sharper and sharper to keep pointing at the target. In the world of automation, turning means wasting (due to the time and energy efficiency). Instead of pointing at the target, the pursuer points its velocity vector ahead of the target, aiming at a predicted future position where their paths will intersect.

\[ a = N \cdot V_c \cdot \dot{\lambda} \]

Where:

\( a \) is the commanded turning acceleration - angular velocity

\( N \) is a constant called the "navigation ratio" (usually between 3 and 5) - how sensitive the chase is due to the direction of the target.

\( V_c \) is the closing velocity - the relative velocity between chaser and target.

\( \dot{\lambda} \) Line-of-Sight (LOS) Angular Rate - the angular velocity between the chaser and target. Play the most important role aka the feedback factor.

8.2 Compare:

Pure Pursuit Curve:

Simple strategy matching the heading direction to the LOS. Give the fastest decision for the device to chase. Using only kinematic path tracing, keep the mathematical difficulty as simple as possible. However, this will require the chaser to have a much higher velocity to catch the target dude to the position follower curve.

Lead Pursuit / Proportional Navigation:

The method used the change in direction of the target to predict a more optimal path. This allows the chaser to reduce the chasing path length, no need for a large speed differencerence, and faster chase regulation with an intelligent target. However, it requires a higher mathematical difficulty to predict the target's path.

9. References

[1] P. Bouguer, Sur les courbes du chien (On the pursuit curve), Mémoires de l’Académie Royale des Sciences, Paris, 1732.

[2] P. J. Nahin, Chases and Escapes: The Mathematics of Pursuit and Evasion, Princeton University Press, 2007.

[3] G. F. Simmons, Differential Equations with Applications and Historical Notes, 3rd ed., CRC Press, 2017.

[4] J. C. Barton and C. J. Eliezer, “On pursuit curves,” Journal of the Australian Mathematical Society (Series B), vol. 41, pp. 358–371, 2000.

[5] Z. Malacka, “Pursuit curves and ordinary differential equations,” Communications, University of Žilina, 2012.

[6] M. Lloyd, “Pursuit curves,” Academic Forum, vol. 24, 2006.

[7] G. Boole, A Treatise on Differential Equations, Macmillan and Co., Cambridge, 1859.

[8] W. H. Press et al., Numerical Recipes: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007.