Department of Mathematics

The Chinese University of Hong Kong

HSMMC Math Modelling Advanced Training Workshop

May 16, 2026

A first-order ordinary differential equation only involves the first derivative of an unknown function with respect to its independent variable. These equations are widely used in mathematics, physics, engineering and many other scientific fields. This paper examines the mathematical form and solution of two classic applications of linear first-order differential equations:

  • Population growth and decay
  • Mixing problems

1 Population Growth and Decay

1.1 Exponential Change

In many real-life situations, how fast a quantity \(y(t)\) changes over time \(t\) depends on its current value. Examples include radioactive material losing mass and population size changing. These processes are called exponential change.

We want to find the function \(y(t)\) that solves this initial value problem (IVP): \[\begin{align} \frac{dy}{dt} &= ky & \quad \quad (1) \\ y(0) &= y_0 & \quad \quad (2) \\ \end{align}\] where \(k\) is a constant.

  • If \(y>0\) and the quantity grows, \(k > 0\).
  • If \(y>0\) and the quantity shrinks, \(k < 0\).
  • If \(y_0 = 0\), then \(y(t) = 0\) for all time \(t\) (a simple solution to Equation (1)).

When \(y \neq 0\), we rearrange Equation (1) to separate variables: \[ \frac{1}{y} \frac{dy}{dt} = k. \]

Integrate both sides with respect to \(t\): \[ \int \frac{1}{y} \frac{dy}{dt} \, dt = \int k \, dt. \]

Using \(dy = \frac{dy}{dt} dt\) on the left side: \[ \int \frac{1}{y} \, dy = \int k \, dt. \]

Calculating the integrals gives: \[ \ln|y| = kt + C, \] where \(C\) is an integration constant. We use exponents to solve for \(y\): \[\begin{align*} |y| &= e^{kt + C} \\ y &= \pm e^{kt} e^{C}. \end{align*}\]

Let \(A = \pm e^{C}\). The general solution is: \[ y(t) = A e^{kt}. \] If we set \(A=0\), this formula also includes the simple solution \(y(t)=0\).

Now we find \(A\) using the initial condition (2). At \(t=0\): \[ y_0 = A e^{0} = A. \]

Substitute \(A = y_0\) to get the final solution: \[ y(t) = y_0 e^{kt}. \quad \quad (3) \]

Equation (3) is exponential growth when \(k > 0\) and exponential decay when \(k < 0\). The constant \(k\) is called the rate constant.

A key feature of radioactive decay is half-life \(T\): the time it takes for a substance to become half its starting amount. For decay (\(k<0\)), set \(y(T) = \frac{1}{2}y_0\) in Equation (3): \[ \frac{1}{2}y_0 = y_0 e^{-kT}. \]

Divide both sides by \(y_0\) (\(y_0 \neq 0\)): \[ \frac{1}{2} = e^{-kT}. \]

Take the natural logarithm of both sides: \[\begin{align*} \ln\left(\frac{1}{2}\right) &= -kT \\ -\ln 2 &= -kT. \end{align*}\]

Solve for \(T\): \[ T = \frac{\ln 2}{k}. \quad \quad (4) \]

1.1.1 Forward Euler Formula

The forward Euler method (forward scheme) is a simple numerical technique to solve the ordinary differential equation (ODE): \[ \frac{dy}{dt} = ky,\quad y(0) = y_0 \]

Forward Euler Formula \[ y_{n+1} = y_n + dt \cdot k \cdot y_n \] where

  • \(y_n\) = current value
  • \(dt\) = time step
  • \(y_{n+1}\) = next value

This code:

  • Implements the forward Euler scheme
  • Compares numerical results with the exact analytical solution (\(y(t) = y_0 e^{kt}\))
  • Plots both solutions for visualization
  • Works for growth (\(k>0\)) and decay (\(k<0\))
import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# User-defined parameters
# --------------------------
y0 = 100        # Initial value (e.g., population, mass)
k = 0.1         # Rate constant: k>0 = growth, k<0 = decay
t_end = 20      # Total simulation time
dt = 0.5        # Time step (smaller = more accurate)

# --------------------------
# Step 1: Create time array
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# Step 2: Forward Euler (Forward Scheme)
# --------------------------
y_euler = np.zeros(n_steps)
y_euler[0] = y0  # Set initial condition

# Iterate using forward scheme
for i in range(n_steps - 1):
    y_euler[i+1] = y_euler[i] + dt * k * y_euler[i]

# --------------------------
# Step 3: Exact analytical solution
# --------------------------
y_exact = y0 * np.exp(k * t)

# --------------------------
# Step 4: Calculate half-life (for decay only)
# --------------------------
if k < 0:
    half_life = np.log(2) / abs(k)
    print(f"=== Radioactive Decay ===")
    print(f"Half-life T = {half_life:.2f} time units")
else:
    print(f"=== Exponential Growth ===")
## === Exponential Growth ===
# --------------------------
# Step 5: Plot results
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, y_euler, 'o-', label='Forward Euler (Numerical)', linewidth=2, markersize=4)
plt.plot(t, y_exact, '-', label='Exact Solution (Analytical)', linewidth=2)
plt.xlabel('Time t')
plt.ylabel('Quantity y(t)')
plt.title(f'Exponential Change (k={k}, y0={y0})')
plt.legend()
plt.grid(True)
plt.show()

# Print first 5 values for verification
print("\nFirst 5 time points:")
## 
## First 5 time points:
print("t\tEuler\t\tExact")
## t    Euler       Exact
for i in range(5):
    print(f"{t[i]:.1f}\t{y_euler[i]:.2f}\t\t{y_exact[i]:.2f}")
## 0.0  100.00      100.00
## 0.5  105.00      105.13
## 1.0  110.25      110.52
## 1.5  115.76      116.18
## 2.0  121.55      122.14

1.1.2 Backward Euler Formula

The Backward Euler method (implicit/backward scheme) is an implicit numerical method for solving the ODE: \[ \frac{dy}{dt} = ky,\quad y(0) = y_0 \]

Backward Euler Formula \[ y_{n+1} = y_n + dt\cdot k\cdot y_{n+1} \]

We rearrange it to get an explicit update rule (solved for \(y_{n+1}\)): \[ y_{n+1} = \frac{y_n}{1 - k\cdot dt} \]

The backward Euler method is implicit (uses \(y_{n+1}\) on both sides), but for the linear ODE \(y' = ky\), we can solve it directly: \[\begin{align*} y_{n+1} &= y_n + k\,dt\,y_{n+1} \\ y_{n+1} - k\,dt\,y_{n+1} &= y_n \\ y_{n+1}(1 - k\,dt) &= y_n \\ y_{n+1} &= \frac{y_n}{1 - k\,dt} \end{align*}\]

This code:

  • Implements the backward Euler scheme
  • Compares with the exact analytical solution
  • Plots results
  • Works for exponential growth (\(k>0\)) and decay (\(k<0\))
import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# Simulation Parameters
# --------------------------
y0 = 100        # Initial value
k = 0.1        # Rate constant (k>0=growth, k<0=decay)
t_end = 20      # Total simulation time
dt = 0.5        # Time step

# --------------------------
# Time array
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# Backward Euler (Backward Scheme)
# --------------------------
y_backward = np.zeros(n_steps)
y_backward[0] = y0

# Iterate with backward Euler formula
for i in range(n_steps - 1):
    # Rearranged explicit formula for linear ODE dy/dt = ky
    y_backward[i+1] = y_backward[i] / (1 - k * dt)

# --------------------------
# Exact Analytical Solution
# --------------------------
y_exact = y0 * np.exp(k * t)

# --------------------------
# Half-life calculation (for decay)
# --------------------------
if k < 0:
    half_life = np.log(2) / abs(k)
    print("=== Exponential Decay ===")
    print(f"Half-life T = {half_life:.2f} time units")
else:
    print("=== Exponential Growth ===")
## === Exponential Growth ===
# --------------------------
# Plot Results
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, y_backward, 's-', label='Backward Euler (Numerical)', linewidth=2, markersize=5)
plt.plot(t, y_exact, '-', label='Exact Solution', linewidth=2)
plt.xlabel('Time t')
plt.ylabel('Quantity y(t)')
plt.title(f'Backward Euler - Exponential Change (k={k}, y0={y0})')
plt.legend()
plt.grid(True)
plt.show()

# --------------------------
# Print sample values
# --------------------------
print("\nFirst 5 time points comparison:")
## 
## First 5 time points comparison:
print(f"t\tBackward Euler\tExact")
## t    Backward Euler  Exact
for i in range(5):
    print(f"{t[i]:.1f}\t{y_backward[i]:.2f}\t\t{y_exact[i]:.2f}")
## 0.0  100.00      100.00
## 0.5  105.26      105.13
## 1.0  110.80      110.52
## 1.5  116.64      116.18
## 2.0  122.77      122.14

1.1.3 Central Difference Scheme

For the ODE \[ \frac{dy}{dt} = ky, \] the central difference scheme (2nd-order accurate) uses the formula: \[ \frac{y_{n+1} - y_{n-1}}{2dt} = k y_n. \]

Rearranged for the update rule: \[ y_{n+1} = y_{n-1} + 2\cdot dt\cdot k\cdot y_n. \]

This is a multi-step method (needs two initial values: \(y_0\) and \(y_1\)). We use Forward Euler to compute \(y_1\) for initialization.

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# Simulation Parameters
# --------------------------
y0 = 100        # Initial value y(0)
k = 0.1        # Rate constant (k>0=growth, k<0=decay)
t_end = 20      # Total simulation time
dt = 0.5        # Time step

# --------------------------
# Create time array
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# Central Difference Scheme (Leapfrog)
# --------------------------
y_central = np.zeros(n_steps)
y_central[0] = y0  # Initial condition

# Step 1: Compute y1 using Forward Euler (initialization for multi-step method)
y_central[1] = y_central[0] + dt * k * y_central[0]

# Step 2: Central difference update for all remaining steps
for n in range(1, n_steps - 1):
    y_central[n+1] = y_central[n-1] + 2 * dt * k * y_central[n]

# --------------------------
# Exact Analytical Solution
# --------------------------
y_exact = y0 * np.exp(k * t)

# --------------------------
# Half-life (for decay only)
# --------------------------
if k < 0:
    half_life = np.log(2) / abs(k)
    print("=== Exponential Decay ===")
    print(f"Half-life T = {half_life:.2f} time units")
else:
    print("=== Exponential Growth ===")
## === Exponential Growth ===
# --------------------------
# Plot Results
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, y_central, 'D-', label='Central Difference Scheme', linewidth=2, markersize=5)
plt.plot(t, y_exact, '-', label='Exact Solution', linewidth=2)
plt.xlabel('Time t')
plt.ylabel('Quantity y(t)')
plt.title(f'Central Difference - Exponential Change (k={k}, y0={y0})')
plt.legend()
plt.grid(True)
plt.show()

# --------------------------
# Print comparison table
# --------------------------
print("\nFirst 5 time points:")
## 
## First 5 time points:
print(f"t\tCentral\t\tExact")
## t    Central     Exact
for i in range(5):
    print(f"{t[i]:.1f}\t{y_central[i]:.2f}\t\t{y_exact[i]:.2f}")
## 0.0  100.00      100.00
## 0.5  105.00      105.13
## 1.0  110.50      110.52
## 1.5  116.05      116.18
## 2.0  122.11      122.14

1.2 Radioactive Decay

Consider a radioactive substance with a half-life of \(1720\) years. At time \(t=0\), the mass of the substance is \(6\) grams. We aim to determine:

  • The remaining mass after \(860\) years.
  • The time \(t_1\) at which the remaining mass is \(2.5\) grams.

From the initial condition, we have \(t_0 = 0\) and \(y_0 = 6\) grams. Substituting these values into the general solution (3), the mass as a function of time is: \[ y(t) = 6 e^{-kt}. \quad \quad (5) \] Note the negative sign in the exponent, which explicitly indicates decay (\(k>0\)).

We are given the half-life \(T = 1720\) years. Using the half-life formula (4), we solve for the rate constant \(k\): \[\begin{align*} k &= \frac{\ln 2}{T} \\ &= \frac{\ln 2}{1720}. \end{align*}\]

Substituting this value of \(k\) into Equation (5) yields the specific solution for this substance: \[ y(t) = 6 e^{-\frac{t \ln 2}{1720}}. \quad \quad (6) \]

1.2.1 Mass After 860 Years

To find the remaining mass after \(860\) years, we evaluate \(y(t)\) at \(t=860\): \[\begin{align*} y(860) &= 6 e^{-\frac{860 \ln 2}{1720}} \\ &= 6 e^{-\frac{\ln 2}{2}} \\ &= 6 \left( e^{\ln 2} \right)^{-1/2} \\ &= 6 (2)^{-1/2} \\ &= \frac{6}{\sqrt{2}} \\ &= 3\sqrt{2} \\ &\approx 4.24 \text{ grams}. \end{align*}\] Thus, the mass remaining after 860 years is approximately \(4.24\) grams.

1.2.2 Time to Reach 2.5 Grams

To find the time \(t_1\) when \(y(t_1) = 2.5\) grams, we substitute these values into Equation (6): \[ 2.5 = 6 e^{-\frac{t_1 \ln 2}{1720}}. \] Note that \(2.5 = \frac{5}{2}\), so the equation becomes: \[ \frac{5}{2} = 6 e^{-\frac{t_1 \ln 2}{1720}}. \] Dividing both sides by \(6\): \[ \frac{5}{12} = e^{-\frac{t_1 \ln 2}{1720}}. \] Taking the natural logarithm of both sides: \[ \ln\left(\frac{5}{12}\right) = -\frac{t_1 \ln 2}{1720}. \] Using the logarithm property \(\ln(a/b) = -\ln(b/a)\), we rewrite the left-hand side: \[ -\ln\left(\frac{12}{5}\right) = -\frac{t_1 \ln 2}{1720}. \] Canceling the negative signs and solving for \(t_1\): \[\begin{align*} t_1 &= 1720 \cdot \frac{\ln\left(\frac{12}{5}\right)}{\ln 2} \\ &\approx 2172.42 \text{ years}. \end{align*}\] Therefore, the mass of the substance will reduce to \(2.5\) grams after approximately \(2172.42\) years.

1.3 Population Growth

The population of Zhuhai was \(1,108,000\) in 2010 and \(1,138,000\) in 2011. Assuming the population grows exponentially at a constant rate, we aim to predict the population in the year 2026.

We define \(t=0\) as the year 2010. The initial population is \(y_0 = 1,108,000\). Using the general solution for exponential growth (3), the population at time \(t\) is: \[ y(t) = 1108000 e^{kt}. \quad \quad (7) \]

To determine the growth constant \(k\), we use the population data from 2011, which corresponds to \(t=1\): \[ y(1) = 1138000 = 1108000 e^{k \cdot 1}. \] Dividing both sides by \(1,108,000\): \[ e^k = \frac{1138000}{1108000} = \frac{1138}{1108}. \] Taking the natural logarithm of both sides: \[ k = \ln\left(\frac{1138}{1108}\right) \approx 0.026716. \]

Substituting this value of \(k\) back into Equation (7), the population model becomes: \[ y(t) = 1108000 e^{0.026716 t}. \]

To find the population in 2026, we note that the year 2026 corresponds to \(t = 2026 - 2010 = 15\). Evaluating the model at \(t=15\): \[\begin{align*} y(15) &= 1108000 e^{0.026716 \cdot 15} \\ &\approx 1108000 e^{0.40074} \\ &\approx 1108000 \cdot 1.4924 \\ &\approx 1654165.39. \end{align*}\] Rounding to the nearest whole number, the population of Zhuhai in 2026 is expected to be approximately \(1,654,165\).

2 Mixing Problems

2.1 Statement of the Problem

Consider a well-mixed tank initially containing \(V_0\) liters of a solution with an initial amount \(A_0\) grams of a dissolved chemical. A solution containing the same chemical at a concentration \(c_1\) grams/liter flows into the tank at a constant rate \(r_1\) liters/minute. The uniformly mixed solution flows out of the tank at a constant rate \(r_2\) liters/minute. Due to continuous stirring, the concentration of the chemical is assumed to be uniform throughout the tank at all times. Let \(A(t)\) be the amount of chemical (in grams) and \(V(t)\) be the volume of solution (in liters) in the tank at time \(t\). The concentration of the chemical leaving the tank, \(c_2(t)\), is then: \[ c_2(t) = \frac{A(t)}{V(t)}. \quad \quad (7) \]

Mixing Problem.

Mixing Problem.

2.2 Mathematical Formulation

We derive the differential equations governing \(V(t)\) and \(A(t)\) by analyzing their changes over a small time interval \(\Delta t\).

The change in volume, \(\Delta V\), is equal to the volume of inflow minus the volume of outflow: \[ \Delta V = r_1 \Delta t - r_2 \Delta t = (r_1 - r_2)\Delta t. \label{eq:vol-change} \quad \quad (8) \]

The change in the amount of chemical, \(\Delta A\), is equal to the amount flowing in minus the amount flowing out. The amount flowing in is \(c_1 r_1 \Delta t\). The amount flowing out is approximately \(c_2 r_2 \Delta t = \frac{A(t)}{V(t)} r_2 \Delta t\). Thus: \[ \Delta A \approx c_1 r_1 \Delta t - \frac{A(t)}{V(t)} r_2 \Delta t = \left(c_1 r_1 - \frac{A(t)}{V(t)} r_2\right)\Delta t. \label{eq:chem-change} \quad \quad (9) \]

Dividing both (8) and (9) by \(\Delta t\) gives the average rates of change: \[ \frac{\Delta V}{\Delta t} = r_1 - r_2, \quad \frac{\Delta A}{\Delta t} \approx c_1 r_1 - \frac{A(t)}{V(t)} r_2. \]

To find the instantaneous rates of change, we take the limit as \(\Delta t \to 0\), which converts the approximations into equalities: \[\begin{align} \frac{dV}{dt} &= r_1 - r_2 \label{eq:ode-vol} & \quad \quad (10) \\ \frac{dA}{dt} &= c_1 r_1 - \frac{A(t)}{V(t)} r_2. \label{eq:ode-chem} & \quad \quad (11) \end{align}\]

We solve for \(V(t)\) first. Since \(r_1\) and \(r_2\) are constants, Equation (10) is a simple ordinary differential equation (ODE) that can be integrated directly. The general solution is: \[ V(t) = (r_1 - r_2)t + C, \] where \(C\) is the constant of integration. Using the initial condition \(V(0) = V_0\), we find \(C=V_0\). Thus, the volume at any time \(t\) is: \[ V(t) = (r_1 - r_2)t + V_0. \]

Substituting this expression for \(V(t)\) into Equation (11) yields a linear first-order ODE for \(A(t)\): \[ \frac{dA}{dt} = c_1 r_1 - \frac{r_2}{(r_1 - r_2)t + V_0} A(t). \] Rearranging into the standard form of a linear ODE, we get: \[ \frac{dA}{dt} + \frac{r_2}{(r_1 - r_2)t + V_0} A(t) = c_1 r_1. \label{eq:linear-ode-chem} \quad \quad (12) \] This equation can be solved using an integrating factor, subject to the initial condition \(A(0) = A_0\).

2.3 General Solution of the Linear Mixing Problem ODE

We consider the first-order linear ordinary differential equation (ODE) governing the amount \(A(t)\) of a chemical in a well-mixed tank, derived from conservation of mass: \[ \frac{dA}{dt} + \frac{r_2}{(r_1 - r_2)t + V_0} A(t) = c_1 r_1. \label{main_ode} \quad \quad (13) \] This ODE satisfies the standard form of a linear first-order ODE: \[ \frac{dA}{dt} + P(t)A(t) = Q(t), \] where we explicitly define the coefficient functions: \[ P(t) = \frac{r_2}{(r_1 - r_2)t + V_0}, \quad Q(t) = c_1 r_1. \] The problem is subject to the initial condition: \[ A(0) = A_0. \label{initial_condition} \quad \quad (14) \] All parameters are constant real numbers:

  • \(r_1\) = volumetric inflow rate (volume per unit time),
  • \(r_2\) = volumetric outflow rate (volume per unit time),
  • \(V_0\) = initial volume of solution in the tank at \(t=0\),
  • \(c_1\) = inflow concentration of the chemical (mass per unit volume),
  • \(A_0\) = initial mass of the chemical in the tank at \(t=0\).

2.3.1 Step 1: Compute the Integrating Factor \(\mu(t)\)

The integrating factor for a linear ODE is defined as: \[ \mu(t) = \exp\left( \int P(t) \, dt \right). \] Substitute \(P(t)\) into the formula: \[ \mu(t) = \exp\left( \int \frac{r_2}{(r_1 - r_2)t + V_0} dt \right). \] We evaluate the integral using u-substitution: Let \[ u = (r_1 - r_2)t + V_0 \implies du = (r_1 - r_2)dt \implies dt = \frac{du}{r_1 - r_2}. \] Substitute into the integral: \[ \int \frac{r_2}{u} \cdot \frac{du}{r_1 - r_2} = \frac{r_2}{r_1 - r_2} \int \frac{1}{u} du. \] The antiderivative of \(1/u\) is the natural logarithm: \[ \frac{r_2}{r_1 - r_2} \ln|u| + C = \frac{r_2}{r_1 - r_2} \ln\left[(r_1 - r_2)t + V_0\right]. \] We drop the absolute value because volume is always positive.

Exponentiate to obtain the integrating factor: \[ \mu(t) = \exp\left( \frac{r_2}{r_1 - r_2} \ln\left[(r_1 - r_2)t + V_0\right] \right). \] Using the logarithm identity \(a\ln b = \ln b^a\) and \(\exp(\ln x)=x\): \[ \mu(t) = \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}}. \quad \quad (15) \label{integrating_factor} \]

2.3.2 Step 2: Multiply the ODE by the Integrating Factor

Multiply both sides of Equation (13) by \(\mu(t)\): \[ \mu(t)\frac{dA}{dt} + \mu(t)P(t)A(t) = \mu(t)Q(t). \] By the design of the integrating factor, the left-hand side is the product rule derivative of \(\mu(t)A(t)\): \[ \frac{d}{dt}\left[ \mu(t) A(t) \right] = c_1 r_1 \cdot \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}}. \label{product_form} \quad \quad (16) \]

2.3.3 Step 3: Integrate Both Sides

Integrate Equation (16) with respect to time \(t\): \[ \mu(t) A(t) = c_1 r_1 \int \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}} dt + C, \] where \(C\) is the constant of integration.

Evaluate the integral with the same substitution \(u=(r_1 - r_2)t + V_0\), \(du=(r_1 - r_2)dt\): \[ \int u^{\frac{r_2}{r_1 - r_2}} \cdot \frac{du}{r_1 - r_2} = \frac{1}{r_1 - r_2} \cdot \frac{u^{\frac{r_2}{r_1 - r_2}+1}}{\frac{r_2}{r_1 - r_2}+1} + C. \] Simplify the exponent: \[ \frac{r_2}{r_1 - r_2} + 1 = \frac{r_2 + (r_1 - r_2)}{r_1 - r_2} = \frac{r_1}{r_1 - r_2}. \] Thus the integral simplifies to: \[ \frac{1}{r_1 - r_2} \cdot \frac{u^{\frac{r_1}{r_1 - r_2}}}{\frac{r_1}{r_1 - r_2}} = \frac{1}{r_1} u^{\frac{r_1}{r_1 - r_2}}. \] Substitute back \(u\) and insert into the equation: \[ \mu(t) A(t) = c_1 r_1 \cdot \frac{1}{r_1}\left[(r_1 - r_2)t + V_0\right]^{\frac{r_1}{r_1 - r_2}} + C. \] Cancel \(r_1\): \[ \mu(t) A(t) = c_1 \left[(r_1 - r_2)t + V_0\right]^{\frac{r_1}{r_1 - r_2}} + C. \quad \quad (17) \label{integrated_eq} \]

2.3.4 Step 4: Solve for \(A(t)\) (General Solution)

Substitute \(\mu(t)\) from Equation (15) into (17): \[ \left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}} A(t) = c_1 \left[(r_1 - r_2)t + V_0\right]^{\frac{r_1}{r_1 - r_2}} + C. \] Divide both sides by \(\left[(r_1 - r_2)t + V_0\right]^{\frac{r_2}{r_1 - r_2}}\): \[ A(t) = c_1 \left[(r_1 - r_2)t + V_0\right] + C \cdot \left[(r_1 - r_2)t + V_0\right]^{-\frac{r_2}{r_1 - r_2}}. \label{general_solution} \quad \quad (18) \] This is the general solution of the ODE.

2.3.5 Step 5: Apply the Initial Condition \(A(0)=A_0\)

Evaluate Equation (18) at \(t=0\): \[ A(0) = c_1 V_0 + C \cdot V_0^{-\frac{r_2}{r_1 - r_2}} = A_0. \] Solve for the constant \(C\): \[ C \cdot V_0^{-\frac{r_2}{r_1 - r_2}} = A_0 - c_1 V_0, \] \[ C = (A_0 - c_1 V_0) \cdot V_0^{\frac{r_2}{r_1 - r_2}}. \label{constant_c} \quad \quad (19) \]

2.3.6 Step 6: Final Particular Solution

Substitute \(C\) from Equation (19) into the general solution (18): \[ \begin{aligned} A(t) &= c_1\left[(r_1 - r_2)t + V_0\right] \\ &\quad + (A_0 - c_1 V_0) V_0^{\frac{r_2}{r_1 - r_2}} \left[(r_1 - r_2)t + V_0\right]^{-\frac{r_2}{r_1 - r_2}}. \end{aligned} \] This is the unique explicit solution to the initial value problem, valid for all \(t\) where the volume remains positive.

2.3.7 A Salt Water Mixing Problem

Consider a tub initially containing \(100\) gallons of salt water in which \(50\) pounds of salt are dissolved. A salt water solution containing \(2\) lb/gal of salt flows into the tub at a rate of \(5\) gal/min. The mixture is kept uniform by stirring and flows out of the tub at a rate of \(4\) gal/min. We aim to find:

  • The rate at which salt enters the tub.
  • The volume of salt water in the tub at time \(t\).
  • The rate at which salt leaves the tub.
  • The solution to the initial value problem describing the mixing process.
  • The concentration of salt in the tub after \(25\) minutes.

We identify the given parameters:

  • Initial volume: \(V_0 = 100\) gal
  • Initial amount of salt: \(A_0 = 50\) lb
  • Inflow concentration: \(c_1 = 2\) lb/gal
  • Inflow rate: \(r_1 = 5\) gal/min
  • Outflow rate: \(r_2 = 4\) gal/min

Rate of Salt Entering

The rate at which salt enters the tub is the product of the inflow concentration and the inflow rate: \[ \text{Rate in} = c_1 r_1 = (2 \text{ lb/gal}) \cdot (5 \text{ gal/min}) = 10 \text{ lb/min}. \]

Volume of Solution at Time \(t\)

The volume of solution in the tub changes at a constant rate. Using the formula derived earlier: \[\begin{align*} V(t) &= V_0 + (r_1 - r_2)t \\ &= 100 + (5 - 4)t \\ &= 100 + t \text{ gal}. \end{align*}\]

Rate of Salt Leaving The concentration of salt in the tub at any time \(t\) is \(\frac{A(t)}{V(t)} = \frac{A(t)}{100 + t}\) lb/gal. The rate at which salt leaves the tub is the product of this concentration and the outflow rate: \[ \text{Rate out} = \left( \frac{A(t)}{100 + t} \text{ lb/gal} \right) \cdot (4 \text{ gal/min}) = \frac{4A(t)}{100 + t} \text{ lb/min}. \]

Differential Equation Model

The rate of change of the amount of salt, \(\frac{dA}{dt}\), is the rate in minus the rate out: \[ \frac{dA}{dt} = 10 - \frac{4A}{100 + t}. \] Rearranging into the standard linear ODE form: \[ \frac{dA}{dt} + \frac{4}{100 + t}A = 10. \] Here, \(P(t) = \frac{4}{100 + t}\) and \(Q(t) = 10\).

Solving the Linear ODE

We solve this equation using an integrating factor, \(\mu(t)\): \[\begin{align*} \mu(t) &= e^{\int P(t) dt} \\ &= e^{\int \frac{4}{100 + t} dt} \\ &= e^{4 \ln|100 + t|} \\ &= e^{\ln(100 + t)^4} \\ &= (100 + t)^4. \end{align*}\]

Multiplying both sides of the standard form ODE by \(\mu(t)\): \[ (100 + t)^4 \frac{dA}{dt} + (100 + t)^4 \frac{4}{100 + t}A = 10(100 + t)^4. \] The left-hand side is the derivative of the product \(\mu(t)A(t)\): \[ \frac{d}{dt} \left[ (100 + t)^4 A(t) \right] = 10(100 + t)^4. \] Integrating both sides with respect to \(t\): \[ (100 + t)^4 A(t) = \int 10(100 + t)^4 dt. \] Evaluating the integral: \[ (100 + t)^4 A(t) = 10 \cdot \frac{(100 + t)^5}{5} + C = 2(100 + t)^5 + C. \] Solving for \(A(t)\) gives the general solution: \[ A(t) = 2(100 + t) + \frac{C}{(100 + t)^4}. \]

We apply the initial condition \(A(0) = 50\) lb to find the constant \(C\): \[\begin{align*} 50 &= 2(100 + 0) + \frac{C}{(100 + 0)^4} \\ 50 &= 200 + \frac{C}{100^4} \\ C &= (50 - 200) \cdot 100^4 = -150 \cdot 100^4. \end{align*}\] Substituting \(C\) back into the general solution yields the particular solution: \[ A(t) = 2(100 + t) - \frac{150 \cdot 100^4}{(100 + t)^4}. \]

Concentration After 25 Minutes

First, we find the amount of salt in the tub at \(t=25\) minutes: \[\begin{align*} A(25) &= 2(100 + 25) - \frac{150 \cdot 100^4}{(100 + 25)^4} \\ &= 250 - \frac{150 \cdot 100^4}{125^4} \\ &\approx 250 - 61.44 \\ &\approx 188.56 \text{ lb}. \end{align*}\] The volume of the solution at \(t=25\) minutes is: \[ V(25) = 100 + 25 = 125 \text{ gal}. \] The concentration of salt at this time is: \[ \text{Concentration} = \frac{A(25)}{V(25)} \approx \frac{188.56}{125} \approx 1.51 \text{ lb/gal}. \]

Central Difference

This code implements the central difference (leapfrog) method for the mixing tank ODE: \[ \frac{dA}{dt} = 10 - \frac{4A}{100+t} \] with initial condition \(A(0)=50\).

Central difference is 2nd-order accurate and far more precise than Euler methods. We use forward Euler to initialize the first step (required for multi-step method).

Central Difference Formula Used \[ A_{i+1} = A_{i-1} + 2\Delta t\cdot \frac{dA}{dt}\left(t_i,A_i\right) \]

  • Multi-step method: needs \(A_0\) and \(A_1\)
  • \(A_1\) computed with forward Euler
  • 2nd-order accuracy \(\to\) very close to exact solution

Example Parameters

  • \(A_0 = 50\) lb
  • \(V_0 = 100\) gal
  • \(r_1 = 5,\ r_2 = 4\) gal/min
  • \(c_1 = 2\) lb/gal
  • Solve up to \(t=25\) min
import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# PARAMETERS (MATCH YOUR PROBLEM)
# --------------------------
A0 = 50          # Initial salt amount (lb)
V0 = 100         # Initial volume (gal)
r1 = 5           # Inflow rate (gal/min)
r2 = 4           # Outflow rate (gal/min)
c1 = 2           # Inlet concentration (lb/gal)
t_end = 25       # Simulate up to 25 minutes
dt = 0.1         # Time step (small = accurate)

# --------------------------
# Time array
# --------------------------
t = np.arange(0, t_end + dt, dt)
n = len(t)

# --------------------------
# Volume function V(t) (analytical)
# --------------------------
def volume(t):
    return V0 + (r1 - r2) * t

# --------------------------
# ODE Function: dA/dt = 10 - (4*A)/(100 + t)
# --------------------------
def dA_dt(t_val, A_val):
    V = volume(t_val)
    return c1*r1 - (r2 * A_val) / V

# --------------------------
# CENTRAL DIFFERENCE (LEAPFROG) SCHEME
# --------------------------
A_central = np.zeros(n)
A_central[0] = A0  # Initial condition

# Step 1: Compute A1 with FORWARD EULER (initialize multi-step method)
A_central[1] = A_central[0] + dt * dA_dt(t[0], A_central[0])

# Step 2: Central difference for all remaining steps
for i in range(1, n-1):
    # Leapfrog formula: A_{i+1} = A_{i-1} + 2*dt*dA/dt(t_i, A_i)
    A_central[i+1] = A_central[i-1] + 2 * dt * dA_dt(t[i], A_central[i])

# --------------------------
# EXACT ANALYTICAL SOLUTION (for comparison)
# --------------------------
def exact_solution(t_val):
    V = volume(t_val)
    C = -150 * (100**4)
    return 2 * V + C / (V**4)

A_exact = exact_solution(t)

# --------------------------
# RESULTS AT t=25 min
# --------------------------
t_target = 25
idx = np.where(np.isclose(t, t_target))[0][0]

V25 = volume(t_target)
A25_central = A_central[idx]
A25_exact = exact_solution(t_target)
conc_central = A25_central / V25
conc_exact = A25_exact / V25

# --------------------------
# PRINT OUTPUT
# --------------------------
print("="*60)
## ============================================================
print("MIXING PROBLEM - CENTRAL DIFFERENCE vs EXACT SOLUTION")
## MIXING PROBLEM - CENTRAL DIFFERENCE vs EXACT SOLUTION
print("="*60)
## ============================================================
print(f"Time = {t_target} min")
## Time = 25 min
print(f"Volume V(t) = {V25:.0f} gal")
## Volume V(t) = 125 gal
print(f"\nSalt Amount A(t):")
## 
## Salt Amount A(t):
print(f"Central Difference: {A25_central:.2f} lb")
## Central Difference: 188.56 lb
print(f"Exact Solution:     {A25_exact:.2f} lb")
## Exact Solution:     188.56 lb
print(f"\nConcentration:")
## 
## Concentration:
print(f"Central Difference: {conc_central:.3f} lb/gal")
## Central Difference: 1.508 lb/gal
print(f"Exact Solution:     {conc_exact:.3f} lb/gal")
## Exact Solution:     1.508 lb/gal
print("="*60)
## ============================================================
# --------------------------
# PLOT
# --------------------------
plt.figure(figsize=(10,5))
plt.plot(t, A_central, 'o-', label='Central Difference', markersize=2, linewidth=1.5)
plt.plot(t, A_exact, '-', label='Exact Solution', linewidth=2)
plt.xlabel('Time t (minutes)')
plt.ylabel('Amount of Salt A(t) (lb)')
plt.title('Mixing Tank Problem: Central Difference Scheme')
plt.legend()
plt.grid(True)
plt.show()

2.3.8 Example: Fertilizer Mixing Problem with a Maximum

A tub initially contains \(100\) gallons of pure water. A solution containing \(1\) lb/gal of fertilizer flows into the tub at a rate of \(1\) gal/min. The mixture is kept uniform by stirring and is pumped out of the tub at a rate of \(3\) gal/min. We aim to find the maximum amount of fertilizer in the tub and the time required to reach this maximum.

Let \(A(t)\) be the amount (in pounds) of fertilizer in the tub at time \(t\). We know that \(A(0) = 0\).

Volume of Solution at Time \(t\)

The volume of solution in the tub changes at a constant rate: \[\begin{align*} V(t) &= V_0 + (r_1 - r_2)t \\ &= 100 + (1 - 3)t \\ &= 100 - 2t \text{ gal}. \end{align*}\] The tub will be empty when \(V(t)=0\), which occurs at \(t=50\) minutes. Thus, the model is valid for \(0 \le t < 50\).

Differential Equation Model

The rate at which fertilizer enters the tub is: \[ \text{Rate in} = c_1 r_1 = (1 \text{ lb/gal}) \cdot (1 \text{ gal/min}) = 1 \text{ lb/min}. \] The rate at which fertilizer leaves the tub is: \[ \text{Rate out} = \left( \frac{A(t)}{V(t)} \text{ lb/gal} \right) \cdot (3 \text{ gal/min}) = \frac{3A(t)}{100 - 2t} \text{ lb/min}. \] The rate of change of the amount of fertilizer, \(\frac{dA}{dt}\), is the rate in minus the rate out: \[ \frac{dA}{dt} = 1 - \frac{3A}{100 - 2t}. \] Rearranging into the standard linear ODE form: \[ \frac{dA}{dt} + \frac{3}{100 - 2t}A = 1. \] Here, \(P(t) = \frac{3}{100 - 2t}\) and \(Q(t) = 1\).

Solving the Linear ODE

We solve this equation using an integrating factor, \(\mu(t)\): \[\begin{align*} \mu(t) &= e^{\int P(t) dt} \\ &= e^{\int \frac{3}{100 - 2t} dt} \\ &= e^{-\frac{3}{2} \ln|100 - 2t|} \\ &= e^{\ln(100 - 2t)^{-3/2}} \\ &= (100 - 2t)^{-3/2}. \end{align*}\]

Multiplying both sides of the standard form ODE by \(\mu(t)\): \[ (100 - 2t)^{-3/2} \frac{dA}{dt} + (100 - 2t)^{-3/2} \frac{3}{100 - 2t}A = (100 - 2t)^{-3/2}. \] The left-hand side is the derivative of the product \(\mu(t)A(t)\): \[ \frac{d}{dt} \left[ (100 - 2t)^{-3/2} A(t) \right] = (100 - 2t)^{-3/2}. \] Integrating both sides with respect to \(t\): \[ (100 - 2t)^{-3/2} A(t) = \int (100 - 2t)^{-3/2} dt. \] Evaluating the integral using substitution \(u=100-2t\), \(du=-2dt\): \[\begin{align*} \int (100 - 2t)^{-3/2} dt &= -\frac{1}{2} \int u^{-3/2} du \\ &= -\frac{1}{2} \left( \frac{u^{-1/2}}{-1/2} \right) + C \\ &= u^{-1/2} + C \\ &= (100 - 2t)^{-1/2} + C. \end{align*}\] Thus, we have: \[ (100 - 2t)^{-3/2} A(t) = (100 - 2t)^{-1/2} + C. \] Solving for \(A(t)\) gives the general solution: \[ A(t) = (100 - 2t) + C(100 - 2t)^{3/2}. \]

We apply the initial condition \(A(0) = 0\) to find the constant \(C\): \[\begin{align*} 0 &= (100 - 0) + C(100 - 0)^{3/2} \\ 0 &= 100 + C(100^{3/2}) \\ 0 &= 100 + C(1000) \\ C &= -\frac{100}{1000} = -\frac{1}{10}. \end{align*}\] Substituting \(C\) back into the general solution yields the particular solution: \[ A(t) = (100 - 2t) - \frac{1}{10}(100 - 2t)^{3/2}. \]

Finding the Maximum

To find the maximum amount of fertilizer, we set the derivative \(\frac{dA}{dt}\) to zero. First, we compute the derivative of our solution: \[\begin{align*} \frac{dA}{dt} &= \frac{d}{dt} \left[ (100 - 2t) - \frac{1}{10}(100 - 2t)^{3/2} \right] \\ &= -2 - \frac{1}{10} \cdot \frac{3}{2}(100 - 2t)^{1/2}(-2) \\ &= -2 + \frac{3}{10}\sqrt{100 - 2t}. \end{align*}\]

Setting \(\frac{dA}{dt} = 0\): \[\begin{align*} -2 + \frac{3}{10}\sqrt{100 - 2t} &= 0 \\ \frac{3}{10}\sqrt{100 - 2t} &= 2 \\ \sqrt{100 - 2t} &= \frac{20}{3} \\ 100 - 2t &= \left(\frac{20}{3}\right)^2 = \frac{400}{9} \\ 2t &= 100 - \frac{400}{9} = \frac{900 - 400}{9} = \frac{500}{9} \\ t &= \frac{250}{9} \approx 27.78 \text{ minutes}. \end{align*}\]

To find the maximum amount of fertilizer, we evaluate \(A(t)\) at \(t = \frac{250}{9}\): \[\begin{align*} A\left(\frac{250}{9}\right) &= \left(100 - 2\left(\frac{250}{9}\right)\right) - \frac{1}{10}\left(100 - 2\left(\frac{250}{9}\right)\right)^{3/2} \\ &= \left(100 - \frac{500}{9}\right) - \frac{1}{10}\left(100 - \frac{500}{9}\right)^{3/2} \\ &= \left(\frac{400}{9}\right) - \frac{1}{10}\left(\frac{400}{9}\right)^{3/2} \\ &= \frac{400}{9} - \frac{1}{10} \left(\frac{20}{3}\right)^3 \\ &= \frac{400}{9} - \frac{1}{10} \left(\frac{8000}{27}\right) \\ &\approx 44.44 - 29.63 \\ &\approx 14.81 \text{ lb}. \end{align*}\] Thus, the maximum amount of fertilizer in the tub is approximately \(14.8\) lb, and it occurs at approximately \(27.78\) minutes.

Central Difference

This code implements the central difference (leapfrog) method to solve your fertilizer mixing ODE: \[ \frac{dA}{dt} = 1 - \frac{3A}{100-2t}, \quad A(0)=0. \]

It validates your analytical solution:

  • Max fertilizer \(\approx 14.81\) lb
  • Time at max \(\approx 27.78\) min

Central difference is 2nd-order accurate and requires forward Euler to initialize the first step.

Central Difference Update Rule: \[ A_{i+1} = A_{i-1} + 2\Delta t \cdot \frac{dA}{dt}\left(t_i, A_i\right). \]

This is a multi-step method - we use forward Euler to compute \(A_1\) before starting the leapfrog iteration.

The code perfectly matches your worked example and verifies the maximum amount and time.

import numpy as np
import matplotlib.pyplot as plt

# --------------------------
# Problem Parameters
# --------------------------
A0 = 0            # Initial fertilizer (lb)
V0 = 100          # Initial volume (gal)
r1 = 1            # Inflow rate (gal/min)
r2 = 3            # Outflow rate (gal/min)
c1 = 1            # Inlet concentration (lb/gal)
t_end = 49.9      # Avoid t=50 (empty tank)
dt = 0.01         # Small time step for accuracy

# --------------------------
# Time array
# --------------------------
t = np.arange(0, t_end + dt, dt)
n_steps = len(t)

# --------------------------
# Volume function V(t)
# --------------------------
def volume(t_val):
    return V0 + (r1 - r2) * t_val

# --------------------------
# ODE: dA/dt = 1 - (3A)/(100 - 2t)
# --------------------------
def dA_dt(t_val, A_val):
    V = volume(t_val)
    return 1 - (3 * A_val) / V

# --------------------------
# CENTRAL DIFFERENCE (LEAPFROG) SCHEME
# --------------------------
A_central = np.zeros(n_steps)
A_central[0] = A0

# Initialize first step with Forward Euler (required for multi-step)
A_central[1] = A_central[0] + dt * dA_dt(t[0], A_central[0])

# Central difference loop
for i in range(1, n_steps - 1):
    A_central[i+1] = A_central[i-1] + 2 * dt * dA_dt(t[i], A_central[i])

# --------------------------
# EXACT ANALYTICAL SOLUTION
# --------------------------
def exact_A(t_val):
    V = volume(t_val)
    return V - (1/10) * (V ** 1.5)

A_exact = exact_A(t)

# --------------------------
# Find Maximum (from numerical solution)
# --------------------------
max_A = np.max(A_central)
t_max = t[np.argmax(A_central)]

# --------------------------
# Print Results
# --------------------------
print("=" * 70)
## ======================================================================
print("FERTILIZER MIXING PROBLEM - CENTRAL DIFFERENCE")
## FERTILIZER MIXING PROBLEM - CENTRAL DIFFERENCE
print("=" * 70)
## ======================================================================
print(f"Maximum fertilizer amount: {max_A:.2f} lb")
## Maximum fertilizer amount: 14.81 lb
print(f"Time at maximum:          {t_max:.2f} min")
## Time at maximum:          27.77 min
print("\nCompare to analytical:")
## 
## Compare to analytical:
print(f"Analytical max:    14.81 lb")
## Analytical max:    14.81 lb
print(f"Analytical time:   27.78 min")
## Analytical time:   27.78 min
print("=" * 70)
## ======================================================================
# --------------------------
# Plot
# --------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, A_central, '.-', label='Central Difference', markersize=1)
plt.plot(t, A_exact, label='Exact Solution', linewidth=2)
plt.axvline(t_max, color='red', linestyle='--', label=f'Max at t={t_max:.2f} min')
plt.xlabel('Time t (min)')
plt.ylabel('Amount of Fertilizer A(t) (lb)')
plt.title('Central Difference Scheme - Fertilizer Mixing Problem')
plt.legend()
plt.grid(True)
plt.show()