Introduction
The integrate.Radau class in SciPy is a powerful method for solving ordinary differential equations (ODEs) with superior accuracy for stiff problems. This tutorial takes you through the essentials of using the Radau method with various examples, starting from basic implementations to more advanced cases that typify real-world scenarios. By the end of this guide, you’ll have a comprehensive understanding of how to leverage Radau within SciPy for your numerical integration tasks.
Prerequisites
- Understanding of ordinary differential equations (ODEs).
- Basic Python programming skills.
- SciPy library installed in your Python environment.
Example #1: Basic Exponential Decay
Our first example demonstrates the application of integrate.Radau to a simple exponential decay problem, represented by the differential equation dy/dt = -ky, where k is a constant.
from scipy.integrate import solve_ivp
# Define the differential equation
def exponential_decay(t, y): return -0.5 * y
# Initial condition
y0 = [10]
# Time span
t_span = (0, 10)
# Solve the equation
sol = solve_ivp(exponential_decay, t_span, y0, method='Radau')
# Output
print('Solution at t=10:', sol.y[0])
Output:
Solution at t=10: [10. 9.67115014 6.9224442 3.74838961 2.02969129 1.09904443
0.59511447 0.32224469 0.1744902 0.09448357 0.06738269]This simple usage showcases how to integrate a differential equation over a specified time interval using the Radau method.
Example #2: Adding a Jacobian for Performance
For systems of equations or more complex dynamics, providing the Jacobian matrix can significantly enhance the performance of the integration. The following example illustrates how to include a Jacobian in a model describing a chemical reaction.
from scipy.integrate import solve_ivp
import numpy as np
# Define the system of ODEs
def chemical_reaction(t, y):
return [
-0.04 * y[0] + 1e4 * y[1] * y[2],
0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1]**2,
3e7 * y[1]**2
]
# Jacobian matrix
def jacobian(t, y):
return np.array([
[-0.04, 1e4 * y[2], 1e4 * y[1]],
[0.04, -1e4 * y[2] - 2 * 3e7 * y[1], -1e4 * y[1]],
[0, 2 * 3e7 * y[1], 0]
])
# Initial conditions
y0 = [1, 0, 0]
# Time span
t_span = (0, 100000)
# Solve the system of ODEs
sol = solve_ivp(chemical_reaction, t_span, y0, method='Radau', jac=jacobian)
# Print results
for i, temp in enumerate(sol.t):
print('Time:', temp, 'Y:', sol.y[:, i])
Output:
Time: 0.0 Y: [1. 0. 0.]
Time: 0.0007066875446339434 Y: [9.99971733e-01 2.37137478e-05 4.55318401e-06]
Time: 0.0014721273912960131 Y: [9.99941120e-01 3.36877365e-05 2.51921817e-05]
Time: 0.002617698972446721 Y: [9.99895320e-01 3.62569765e-05 6.84231398e-05]
Time: 0.0038099051153662345 Y: [9.99847678e-01 3.64772143e-05 1.15845034e-04]
Time: 0.007403197788978978 Y: [9.99704224e-01 3.64693124e-05 2.59306773e-04]
Time: 0.019926790101604357 Y: [9.99205871e-01 3.63782129e-05 7.57750591e-04]
Time: 0.09159974354695757 Y: [9.96401273e-01 3.58634001e-05 3.56286323e-03]
Time: 0.31587606451877004 Y: [9.88109520e-01 3.43776903e-05 1.18561018e-02]
Time: 0.9639800232767326 Y: [9.67486558e-01 3.09104303e-05 3.24825320e-02]
Time: 1.612083982034695 Y: [9.50502854e-01 2.82988392e-05 4.94688470e-02]
Time: 2.9557170244645947 Y: [9.22677906e-01 2.44837763e-05 7.72976104e-02]
Time: 4.7281478552432326 Y: [8.95131103e-01 2.12411103e-05 1.04847656e-01]
Time: 7.4053045388200704 Y: [8.64275940e-01 1.81735835e-05 1.35705886e-01]
Time: 11.04974163954623 Y: [8.33393538e-01 1.56170330e-05 1.66590845e-01]
Time: 16.622814534799943 Y: [7.98937779e-01 1.32681771e-05 2.01048953e-01]
Time: 24.351444486358574 Y: [7.64261750e-01 1.13318870e-05 2.35726918e-01]
Time: 36.29437078896136 Y: [7.25604432e-01 9.57535434e-06 2.74385992e-01]
Time: 53.04010390580073 Y: [6.86701648e-01 8.13214765e-06 3.13290219e-01]
Time: 79.0939628005826 Y: [6.43544099e-01 6.82931918e-06 3.56449072e-01]
Time: 115.74535358063251 Y: [6.00514720e-01 5.76213752e-06 3.99479518e-01]
Time: 172.75141656820014 Y: [5.53453899e-01 4.80280090e-06 4.46541299e-01]
Time: 252.48292981773693 Y: [5.07491649e-01 4.02225999e-06 4.92504329e-01]
Time: 375.561902270049 Y: [4.58380529e-01 3.32397220e-06 5.41616147e-01]
Time: 545.1954978375046 Y: [4.11885329e-01 2.76201914e-06 5.88111909e-01]
Time: 803.7755419584241 Y: [3.63653877e-01 2.26168037e-06 6.36343862e-01]
Time: 1153.7900233906576 Y: [3.19589735e-01 1.86331096e-06 6.80408402e-01]
Time: 1677.2864318807187 Y: [2.75530353e-01 1.51175643e-06 7.24468135e-01]
Time: 2371.2090026194155 Y: [2.36798010e-01 1.23502498e-06 7.63200755e-01]
Time: 3387.645180895694 Y: [1.99516067e-01 9.93251503e-07 8.00482939e-01]
Time: 4706.760656280608 Y: [1.67956525e-01 8.05095567e-07 8.32042670e-01]
Time: 6600.997983353201 Y: [1.38641590e-01 6.42379952e-07 8.61357767e-01]
Time: 9012.168847499417 Y: [1.14660060e-01 5.17135093e-07 8.85339423e-01]
Time: 12416.139522510943 Y: [9.30559025e-02 4.09856520e-07 9.06943688e-01]
Time: 16679.56360295253 Y: [7.58804693e-02 3.28097986e-07 9.24119203e-01]
Time: 22618.232526381784 Y: [6.07811590e-02 2.58644089e-07 9.39218582e-01]
Time: 29964.9488345872 Y: [4.90406811e-02 2.06146450e-07 9.50959113e-01]
Time: 40099.20600008357 Y: [3.89060664e-02 1.61842329e-07 9.61093772e-01]
Time: 52528.53993897204 Y: [3.11521621e-02 1.28564922e-07 9.68847709e-01]
Time: 69564.10149634449 Y: [2.45438494e-02 1.00618148e-07 9.75456050e-01]
Time: 92064.83952727004 Y: [1.92228659e-02 7.83801007e-08 9.80777056e-01]
Time: 100000.0 Y: [1.78655726e-02 7.27460704e-08 9.82134355e-01]In this example, the precise definition of the Jacobian matrix improves the integration speed and accuracy.
Example #3: Handling Events in Integration
The Radau solver can also effectively handle events during integration. Events are conditions that, when met, can trigger actions such as halting the integration. This is illustrated in the following example, where the integration is stopped at a predefined value of the dependent variable.
from scipy.integrate import solve_ivp
# Event function
def event(t, y): return y[0] - 0.1
event.terminal = True
event.direction = -1
# Differential equation
def model(t, y): return -0.5 * y
# Initial condition and time span
y0 = [2]
t_span = (0, 5)
# Solve with event
sol = solve_ivp(model, t_span, y0, method='Radau', events=event)
# Output
print('Integration stopped at t:', sol.t_events[0])
This capability is especially useful for stopping the integration process once the solution reaches a specific state, without the need for manual monitoring.
Example #4: Solving Systems with Complex Coefficients
The integrate.Radau class can effectively handle systems of ODEs with complex coefficients, common in electromagnetic and quantum mechanical problems. Here, we demonstrate a simple harmonic oscillator with a complex spring constant.
from scipy.integrate import solve_ivp
import numpy as np
# Oscillator equation
def oscillator(t, y): return [y[1], -0.15 - 0.1j * y[0]]
# Initial condition
y0 = [1 + 0j, 0 + 0j]
# Time span
t_span = (0, 50)
# Solve
sol = solve_ivp(oscillator, t_span, y0, method='Radau')
# Show real part of the solution
import matplotlib.pyplot as plt
plt.plot(sol.t, sol.y[0].real)
plt.title('Real part of Oscillator Displacement')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.show()
Handling complex coefficients powers solutions to a wide range of scientific and engineering challenges.
Example #5: Coupled Differential Equations with External Inputs
Lastly, we explore how integrate.Radau can manage coupled differential equations affected by external inputs or controls, a common scenario in control systems and environmental modeling. This example simulates a system influenced by an external temperature change over time.
from scipy.integrate import solve_ivp
import numpy as np
# System of equations with external input
def thermal_system(t, y, T_external):\n dydt = [y[1], -0.5*y[1] + 0.5*(T_external(t) - y[0])]\n return dydt
# External temperature function
def T_external(t): return 20 + 10*np.sin(t)
# Wrap the system for solve_ivp, passing external function as argument
def wrapped_system(t, y): return thermal_system(t, y, T_external)
# Initial conditions
y0 = [20, 0]
# Time span
t_span = (0, 24)
# Integration
sol = solve_ivp(wrapped_system, t_span, y0, method='Radau')
# Result
plt.plot(sol.t, sol.y[0])
plt.title('Temperature Response to External Changes')
plt.xlabel('Time (s)')
plt.ylabel('Temperature (C)')
plt.show()
This example highlights the flexibility of integrate.Radau in solving complex systems with time-dependent external influences.
Conclusion
The integrate.Radau class in SciPy is a versatile tool for solving stiff differential equations, offering high accuracy and the ability to handle complex scenarios. Through the examples provided, it’s clear that whether you’re dealing with simple models or intricate systems impacted by external conditions, the Radau method can be a crucial component in your numerical computation arsenal.
