Introduction

Ordinary Differential Equations (ODEs) are mathematical tools used to model dynamic systems, and solving them numerically is pivotal in various scientific domains. In this technical blog, we delve into the numerical solution of Robertson’s ODE, a system of three coupled ODEs. The Python code provided utilizes both SciPy’s built-in solver and a custom implementation based on Euler’s method with a modified Euler corrector. Let’s explore the underlying theory and significance of Robertson’s ODE and the numerical methods employed.

Robertson’s ODE

The system of ODEs represents a chemical reaction network involving three species ((y_1), (y_2), and (y_3)) with intricate interactions. Mathematically, the ODEs are defined as follows:

\[\frac{dy_1}{dt} = -0.04 \cdot y_1 + 10000 \cdot y_2 \cdot y_3 \\ \frac{dy_2}{dt} = 0.04 \cdot y_1 - 10000 \cdot y_2 \cdot y_3 - 30000000 \cdot y_2^2 \\ \frac{dy_3}{dt} = 30000000 \cdot y_2^2\]

Here, \(y_1\), \(y_2\), and \(y_3\) represent concentrations of different chemical species, and the derivatives \(\frac{dy_i}{dt}\) describe their rates of change.

Numerical Solution using SciPy

The code utilizes the LSODA method from scipy.integrate.solve_ivp for the numerical solution. LSODA is a versatile solver that adapts its step size to ensure accuracy in capturing the dynamics of the system. The resulting solution is then visualized using Matplotlib.

import numpy as np
import matplotlib.pyplot as plt

def robertson(t,y):
  y1 = y[0]
  y2 = y[1]
  y3 = y[2]

  dydt = np.zeros(3)

  dydt[0] = -0.04 * y1 + 10000.0 * y2 * y3
  dydt[1] = 0.04 * y1 - 10000.0 * y2 * y3 - 30000000.0 * y2 * y2
  dydt[2] = 30000000.0 * y2 * y2

  return dydt

from scipy.integrate import solve_ivp
x0=np.asarray([1,0,0])
t0=0
tn=40
tspan = np.array ( [ t0, tn ] )
sol = solve_ivp ( robertson, tspan, x0, method = 'LSODA')
sol

plt.plot ( sol.t, sol.y[0,:], linewidth = 3 )
plt.plot ( sol.t, sol.y[1,:], linewidth = 3 )
plt.plot ( sol.t, sol.y[2,:], linewidth = 3 )
plt.grid ( True )
plt.xlabel ( 't' )
plt.ylabel ( 'y(t)' )
plt.title ( 'Robertson ODE' )
plt.legend ( [ 'y1', 'y2', 'y3' ] );

This approach is convenient for its efficiency and accuracy, especially when dealing with stiff ODEs like those in chemical kinetics.

Numerical Solution using Euler’s Method

While the SciPy solution provides an accurate result, implementing a custom numerical method offers insights into the mechanics of solving ODEs. The Euler method, a simple yet fundamental numerical technique, is employed. The code includes functions for both the predictor and corrector steps, providing a hands-on understanding of the algorithm.

def euler_predictor(f, x, y, h):
    predicted_y = y + h * f(x, y)
    return predicted_y

def modified_euler_corrector(f, x, y, x1, y1_predicted, h):
    epsilon = 0.00000001
    corrected_y = y + 0.5 * h * (f(x, y) + f(x1, y1_predicted))

    while np.sum(np.abs(corrected_y - y1_predicted)) > epsilon:
        y1_predicted = corrected_y
        corrected_y = y + 0.5 * h * (f(x, y) + f(x1, y1_predicted))

    return corrected_y

def solve_differential_equation(f, initial_x, final_x, initial_y, step_size):
    x_values = [initial_x]
    y_values = [initial_y]

    while initial_x < final_x:
        next_x = initial_x + step_size
        predicted_y = euler_predictor(f, initial_x, initial_y, step_size)
        corrected_y = modified_euler_corrector(f, initial_x, initial_y, next_x, predicted_y, step_size)
        initial_x = next_x
        initial_y = corrected_y
        x_values.append(initial_x)
        y_values.append(initial_y)

    
    return np.asarray(x_values), np.asarray(y_values)


x = 0;
y = np.asarray([1,0,0]);
xn = 40;
h = .00003;
xstor,ystor=solve_differential_equation(robertson,x, xn, y, h);

plt.plot(xstor,ystor);
plt.grid ( True )
plt.xlabel ( 't' )
plt.ylabel ( 'y(t)' )
plt.title ( 'Robertson ODE' )
plt.legend ( [ 'y1', 'y2', 'y3' ] );

This approach allows for a step-by-step exploration of the solution, emphasizing the iterative nature of numerical ODE solvers.

Theoretical Insights

Stiffness in ODEs

The system’s stiffness, a measure of how rapidly some components of the solution change compared to others, is a crucial consideration. The stiffness in Robertson’s ODE arises from the vastly different timescales involved in the reactions. SciPy’s LSODA method excels in handling stiff ODEs by dynamically adjusting the step size.

Euler’s Method and Correctors

Euler’s method, while straightforward, may be limited in accuracy for stiff systems. The introduction of a modified Euler corrector enhances accuracy by refining the solution iteratively. This correction step ensures a more reliable approximation of the true solution.

Conclusion

In this blog, we explored the numerical solution of Robertson’s ODE using both SciPy’s powerful solver and a custom Euler method with a modified corrector. Understanding the theoretical underpinnings, including the concept of stiffness and the mechanics of numerical methods, provides a holistic perspective on solving complex ODEs. Both the convenience of a high-level library and the insights gained from a manual approach contribute to a comprehensive understanding of numerical integration in the context of dynamic systems.