Integration and Differential Equations
SciPy provides powerful numerical integration and differential equation solving capabilities primarily through the scipy.integrate module. This module includes various numerical integration methods, ordinary differential equation (ODE) solvers, partial differential equation (PDE) solving tools, and more. These functions have wide applications in scientific computing, engineering analysis, and physical modeling. This chapter will provide a detailed introduction to using SciPy for numerical integration and differential equation solving.
scipy.integrate Module Overview
The scipy.integrate module includes the following main functions:
- Numerical integration (definite and indefinite integrals)
- Ordinary differential equation solving
- Boundary value problem solving
- Integral equation solving
- Multiple integral calculation
python
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
from scipy.integrate import (
quad, dblquad, tplquad, nquad,
odeint, solve_ivp, solve_bvp,
simpson, trapz, cumtrapz,
fixed_quad, quadrature
)
import warnings
warnings.filterwarnings('ignore')
# Set plot style
plt.style.use('seaborn-v0_8')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# View main functions of integrate module
print("scipy.integrate Main Functions:")
functions = [attr for attr in dir(integrate) if not attr.startswith('_')]
print(f"Total {len(functions)} functions and classes")
print("Integration functions:", [f for f in functions if 'quad' in f or 'int' in f][:10])
print("ODE solvers:", [f for f in functions if 'ode' in f or 'solve' in f][:10])Numerical Integration
1. One-Dimensional Integration
python
# Define test functions
def test_functions():
"""Define various test functions"""
functions = {
'Polynomial': lambda x: x**3 - 2*x**2 + x - 1,
'Trigonometric': lambda x: np.sin(x) * np.cos(x),
'Exponential': lambda x: np.exp(-x**2),
'Rational': lambda x: 1 / (1 + x**2),
'Oscillating': lambda x: np.sin(10*x) / x if x != 0 else 10,
'Discontinuous': lambda x: np.where(x < 0.5, x**2, 2*x - 1)
}
# Analytical solutions (for verification)
analytical_solutions = {
'Polynomial': lambda a, b: (b**4/4 - 2*b**3/3 + b**2/2 - b) - (a**4/4 - 2*a**3/3 + a**2/2 - a),
'Trigonometric': lambda a, b: (np.sin(b)**2 - np.sin(a)**2) / 2,
'Rational': lambda a, b: np.arctan(b) - np.arctan(a)
}
return functions, analytical_solutions
# Get test functions
test_funcs, analytical_sols = test_functions()
# Use different integration methods
integration_methods = {
'quad (Adaptive)': lambda f, a, b: quad(f, a, b),
'fixed_quad (Gaussian)': lambda f, a, b: fixed_quad(f, a, b, n=5),
'quadrature (Adaptive Gaussian)': lambda f, a, b: quadrature(f, a, b)
}
# Integration interval
a, b = 0, 1
print(f"Numerical Integration Comparison (Interval [{a}, {b}]):")
print(f"{'Function':>15} {'Method':>15} {'Numerical Result':>12} {'Error Estimate':>12} {'Analytical':>12} {'Absolute Error':>12}")
print("-" * 90)
for func_name, func in test_funcs.items():
# Calculate analytical solution (if available)
analytical_result = None
if func_name in analytical_sols:
try:
analytical_result = analytical_sols[func_name](a, b)
except:
pass
for method_name, method in integration_methods.items():
try:
result, error = method(func, a, b)
# Calculate error compared to analytical solution
abs_error = ""
if analytical_result is not None:
abs_error = f"{abs(result - analytical_result):.2e}"
analytical_str = f"{analytical_result:.6f}" if analytical_result is not None else "N/A"
print(f"{func_name:>15} {method_name:>15} {result:>12.6f} {error:>12.2e} {analytical_str:>12} {abs_error:>12}")
except Exception as e:
print(f"{func_name:>15} {method_name:>15} {'Failed':>12} {'N/A':>12} {'N/A':>12} {'N/A':>12}")
print()
# Visualize integration process
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
x = np.linspace(0, 1, 1000)
for i, (name, func) in enumerate(test_funcs.items()):
if i >= 6:
break
# Calculate function values
try:
y = [func(xi) for xi in x]
axes[i].plot(x, y, 'b-', linewidth=2, label=f'{name}')
axes[i].fill_between(x, 0, y, alpha=0.3, color='lightblue')
axes[i].set_title(f'{name}')
axes[i].set_xlabel('x')
axes[i].set_ylabel('f(x)')
axes[i].grid(True, alpha=0.3)
axes[i].legend()
# Add integral value annotation
integral_value, _ = quad(func, a, b)
axes[i].text(0.5, max(y)*0.8, f'∫f(x)dx = {integral_value:.4f}',
ha='center', va='center',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
except Exception as e:
axes[i].text(0.5, 0.5, f'Plot Failed\n{str(e)}', ha='center', va='center')
axes[i].set_title(f'{name} (Plot Failed)')
plt.tight_layout()
plt.show()2. Multiple Integration
python
# Double integration examples
def double_integration_examples():
"""Double integration examples"""
# Define bivariate functions
functions_2d = {
'Simple Polynomial': lambda y, x: x*y,
'Gaussian Function': lambda y, x: np.exp(-(x**2 + y**2)),
'Trigonometric': lambda y, x: np.sin(x) * np.cos(y),
'Rational': lambda y, x: 1 / (1 + x**2 + y**2)
}
# Integration regions
regions = {
'Rectangular Region': {'x_bounds': [0, 1], 'y_bounds': [0, 1]},
'Triangular Region': {'x_bounds': [0, 1], 'y_bounds': lambda x: [0, x]},
'Circular Region': {'x_bounds': [-1, 1], 'y_bounds': lambda x: [-np.sqrt(1-x**2), np.sqrt(1-x**2)]}
}
print("Double Integral Calculation:")
print(f"{'Function':>15} {'Region':>15} {'Integral Value':>12} {'Error Estimate':>12}")
print("-" * 60)
results = {}
for func_name, func in functions_2d.items():
results[func_name] = {}
for region_name, region in regions.items():
try:
if callable(region['y_bounds']):
# Variable limits integration
result, error = dblquad(func,
region['x_bounds'][0], region['x_bounds'][1],
lambda x: region['y_bounds'](x)[0],
lambda x: region['y_bounds'](x)[1])
else:
# Constant limits integration
result, error = dblquad(func,
region['x_bounds'][0], region['x_bounds'][1],
region['y_bounds'][0], region['y_bounds'][1])
results[func_name][region_name] = (result, error)
print(f"{func_name:>15} {region_name:>15} {result:>12.6f} {error:>12.2e}")
except Exception as e:
print(f"{func_name:>15} {region_name:>15} {'Failed':>12} {'N/A':>12}")
results[func_name][region_name] = (None, None)
return results
# Execute double integration calculation
double_integral_results = double_integration_examples()
# Triple integral examples
print("\nTriple Integral Calculation:")
# Define ternary function
def triple_function(z, y, x):
return x*y*z
# Calculate triple integral
try:
# Integrate on unit cube
result_3d, error_3d = tplquad(triple_function, 0, 1, 0, 1, 0, 1)
print(f"∭xyz dxdydz (Unit Cube) = {result_3d:.6f} ± {error_3d:.2e}")
# Analytical solution: 1/8
analytical_3d = 1/8
print(f"Analytical Solution = {analytical_3d:.6f}")
print(f"Absolute Error = {abs(result_3d - analytical_3d):.2e}")
except Exception as e:
print(f"Triple integral calculation failed: {e}")
# Visualize double integration
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()
# Create grid
x = np.linspace(-1, 1, 50)
y = np.linspace(-1, 1, 50)
X, Y = np.meshgrid(x, y)
# Visualize different bivariate functions
functions_for_plot = {
'Gaussian Function': lambda x, y: np.exp(-(x**2 + y**2)),
'Trigonometric': lambda x, y: np.sin(x) * np.cos(y),
'Rational': lambda x, y: 1 / (1 + x**2 + y**2),
'Wave Function': lambda x, y: np.sin(np.sqrt(x**2 + y**2))
}
for i, (name, func) in enumerate(functions_for_plot.items()):
Z = func(X, Y)
contour = axes[i].contourf(X, Y, Z, levels=20, cmap='viridis')
axes[i].set_title(f'{name}')
axes[i].set_xlabel('x')
axes[i].set_ylabel('y')
plt.colorbar(contour, ax=axes[i])
# Add integration region markers
circle = plt.Circle((0, 0), 1, fill=False, color='red', linewidth=2, linestyle='--')
axes[i].add_patch(circle)
axes[i].text(0, -1.3, 'Red dashed line: Circular integration region', ha='center', color='red')
plt.tight_layout()
plt.show()3. Comparison of Numerical Integration Methods
python
# Compare precision and efficiency of different numerical integration methods
import time
def compare_integration_methods():
"""Compare different integration methods"""
# Test functions
def smooth_function(x):
return np.exp(-x**2) * np.cos(x)
def oscillatory_function(x):
return np.sin(20*x) / (1 + x**2)
def singular_function(x):
return 1 / np.sqrt(x) if x > 0 else 0
test_functions = {
'Smooth Function': smooth_function,
'Oscillatory Function': oscillatory_function,
'Singular Function': singular_function
}
# Integration methods
methods = {
'quad': lambda f, a, b: quad(f, a, b),
'simpson': lambda f, a, b: (simpson([f(x) for x in np.linspace(a, b, 1001)], dx=(b-a)/1000), 0),
'trapz': lambda f, a, b: (trapz([f(x) for x in np.linspace(a, b, 1001)], dx=(b-a)/1000), 0),
'fixed_quad': lambda f, a, b: fixed_quad(f, a, b, n=10)
}
# Integration intervals
intervals = {
'Smooth Function': (0, 2),
'Oscillatory Function': (0, 1),
'Singular Function': (0.001, 1) # Avoid singularity
}
print("Integration Method Comparison:")
print(f"{'Function':>12} {'Method':>12} {'Result':>12} {'Time(ms)':>10} {'Relative Error':>12}")
print("-" * 70)
for func_name, func in test_functions.items():
a, b = intervals[func_name]
# Use quad as reference
reference_result, _ = quad(func, a, b)
for method_name, method in methods.items():
try:
start_time = time.time()
result, _ = method(func, a, b)
exec_time = (time.time() - start_time) * 1000
relative_error = abs(result - reference_result) / abs(reference_result) if reference_result != 0 else 0
print(f"{func_name:>12} {method_name:>12} {result:>12.6f} {exec_time:>10.2f} {relative_error:>12.2e}")
except Exception as e:
print(f"{func_name:>12} {method_name:>12} {'Failed':>12} {'N/A':>10} {'N/A':>12}")
print()
# Execute method comparison
compare_integration_methods()
# Visualize integration accuracy
def visualize_integration_accuracy():
"""Visualize integration accuracy"""
def test_func(x):
return np.sin(x) * np.exp(-x/2)
# Analytical solution
def analytical_integral(a, b):
def F(x):
return -2/5 * np.exp(-x/2) * (np.sin(x) + 2*np.cos(x))
return F(b) - F(a)
a, b = 0, 2*np.pi
true_value = analytical_integral(a, b)
# Different grid densities
n_points = [5, 10, 20, 50, 100, 200, 500, 1000]
methods_accuracy = {
'Simpson': [],
'Trapezoid': [],
'Quad': []
}
for n in n_points:
x = np.linspace(a, b, n)
y = test_func(x)
dx = (b - a) / (n - 1)
# Simpson method
if n >= 3:
simpson_result = simpson(y, dx=dx)
methods_accuracy['Simpson'].append(abs(simpson_result - true_value))
else:
methods_accuracy['Simpson'].append(np.nan)
# Trapezoid method
trapz_result = trapz(y, dx=dx)
methods_accuracy['Trapezoid'].append(abs(trapz_result - true_value))
# Quad method (reference)
quad_result, _ = quad(test_func, a, b)
methods_accuracy['Quad'].append(abs(quad_result - true_value))
# Plot accuracy comparison
plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
x_plot = np.linspace(a, b, 1000)
y_plot = test_func(x_plot)
plt.plot(x_plot, y_plot, 'b-', linewidth=2, label='f(x) = sin(x)exp(-x/2)')
plt.fill_between(x_plot, 0, y_plot, alpha=0.3)
plt.title('Integrand Function')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(2, 2, 2)
for method, errors in methods_accuracy.items():
if method != 'Quad': # Quad error is too small, display separately
plt.loglog(n_points, errors, 'o-', label=method, linewidth=2, markersize=6)
plt.xlabel('Grid Points')
plt.ylabel('Absolute Error')
plt.title('Integration Accuracy Comparison')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(2, 2, 3)
# Show convergence order
simpson_errors = [e for e in methods_accuracy['Simpson'] if not np.isnan(e)]
trapz_errors = methods_accuracy['Trapezoid']
if len(simpson_errors) > 1:
simpson_order = np.polyfit(np.log(n_points[-len(simpson_errors):]),
np.log(simpson_errors), 1)[0]
plt.text(0.1, 0.8, f'Simpson Convergence Order: {-simpson_order:.2f}', transform=plt.gca().transAxes)
trapz_order = np.polyfit(np.log(n_points), np.log(trapz_errors), 1)[0]
plt.text(0.1, 0.7, f'Trapezoid Convergence Order: {-trapz_order:.2f}', transform=plt.gca().transAxes)
plt.text(0.1, 0.6, f'True Integral Value: {true_value:.6f}', transform=plt.gca().transAxes)
plt.text(0.1, 0.5, f'Quad Result: {quad_result:.6f}', transform=plt.gca().transAxes)
plt.text(0.1, 0.4, f'Quad Error: {abs(quad_result - true_value):.2e}', transform=plt.gca().transAxes)
plt.title('Convergence Analysis')
plt.axis('off')
plt.subplot(2, 2, 4)
# Show final accuracy of different methods
final_errors = {method: errors[-1] for method, errors in methods_accuracy.items()}
methods_list = list(final_errors.keys())
errors_list = list(final_errors.values())
bars = plt.bar(methods_list, errors_list, color=['red', 'green', 'blue'], alpha=0.7)
plt.yscale('log')
plt.ylabel('Absolute Error')
plt.title('Final Accuracy Comparison')
plt.grid(True, alpha=0.3)
# Add numerical labels
for bar, error in zip(bars, errors_list):
plt.text(bar.get_x() + bar.get_width()/2, error, f'{error:.2e}',
ha='center', va='bottom')
plt.tight_layout()
plt.show()
# Execute accuracy visualization
visualize_integration_accuracy()Ordinary Differential Equation Solving
1. Initial Value Problems
python
# Ordinary differential equation initial value problem solving
def ode_examples():
"""Ordinary differential equation examples"""
# Example 1: First-order linear ODE
# dy/dt = -2y + 1, y(0) = 0
# Analytical solution: y(t) = 0.5(1 - exp(-2t))
def linear_ode(t, y):
return -2*y + 1
def linear_analytical(t):
return 0.5 * (1 - np.exp(-2*t))
# Example 2: Nonlinear ODE (Logistic equation)
# dy/dt = ry(1 - y/K), y(0) = y0
def logistic_ode(t, y, r=1, K=10):
return r * y * (1 - y/K)
def logistic_analytical(t, y0=1, r=1, K=10):
return K / (1 + (K/y0 - 1) * np.exp(-r*t))
# Example 3: Oscillator (Simple harmonic motion)
# d²y/dt² + ω²y = 0 => dy/dt = v, dv/dt = -ω²y
def harmonic_oscillator(t, y, omega=2):
y_pos, y_vel = y
return [y_vel, -omega**2 * y_pos]
def harmonic_analytical(t, y0=1, v0=0, omega=2):
A = np.sqrt(y0**2 + (v0/omega)**2)
phi = np.arctan2(-v0/omega, y0)
return A * np.cos(omega*t + phi)
return {
'linear': (linear_ode, linear_analytical, [0], 'dy/dt = -2y + 1'),
'logistic': (logistic_ode, lambda t: logistic_analytical(t), [1], 'dy/dt = ry(1-y/K)'),
'harmonic': (harmonic_oscillator, lambda t: harmonic_analytical(t), [1, 0], 'd²y/dt² + ω²y = 0')
}
# Get ODE examples
ode_examples = ode_examples()
# Solve ODE and compare methods
t_span = (0, 5)
t_eval = np.linspace(0, 5, 100)
print("Ordinary Differential Equation Solving Comparison:")
print(f"{'Equation':>15} {'Method':>15} {'Final Value':>12} {'Analytical':>12} {'Absolute Error':>12} {'Time(ms)':>10}")
print("-" * 85)
# Different solving methods
methods = {
'RK45': 'RK45',
'RK23': 'RK23',
'DOP853': 'DOP853',
'Radau': 'Radau',
'BDF': 'BDF'
}
results = {}
for ode_name, (ode_func, analytical_func, y0, description) in ode_examples.items():
results[ode_name] = {}
# Calculate analytical solution
if ode_name == 'harmonic':
analytical_values = [harmonic_analytical(t) for t in t_eval]
final_analytical = analytical_values[-1]
else:
analytical_values = [analytical_func(t) for t in t_eval]
final_analytical = analytical_values[-1]
for method_name, method in methods.items():
try:
start_time = time.time()
if ode_name == 'logistic':
# ODE with parameters
sol = solve_ivp(lambda t, y: logistic_ode(t, y, r=1, K=10),
t_span, y0, t_eval=t_eval, method=method)
else:
sol = solve_ivp(ode_func, t_span, y0, t_eval=t_eval, method=method)
exec_time = (time.time() - start_time) * 1000
if sol.success:
if ode_name == 'harmonic':
final_numerical = sol.y[0, -1] # Position component
else:
final_numerical = sol.y[0, -1]
abs_error = abs(final_numerical - final_analytical)
results[ode_name][method_name] = {
'solution': sol,
'final_value': final_numerical,
'error': abs_error,
'time': exec_time
}
print(f"{ode_name:>15} {method_name:>15} {final_numerical:>12.6f} {final_analytical:>12.6f} {abs_error:>12.2e} {exec_time:>10.2f}")
else:
print(f"{ode_name:>15} {method_name:>15} {'Failed':>12} {'N/A':>12} {'N/A':>12} {'N/A':>10}")
except Exception as e:
print(f"{ode_name:>15} {method_name:>15} {'Error':>12} {'N/A':>12} {'N/A':>12} {'N/A':>10}")
print()
# Visualize ODE solving results
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
axes = axes.flatten()
plot_idx = 0
for ode_name, (ode_func, analytical_func, y0, description) in ode_examples.items():
if plot_idx >= 3:
break
ax = axes[plot_idx]
# Plot analytical solution
if ode_name == 'harmonic':
analytical_values = [harmonic_analytical(t) for t in t_eval]
ax.plot(t_eval, analytical_values, 'k--', linewidth=3, label='Analytical Solution', alpha=0.8)
else:
analytical_values = [analytical_func(t) for t in t_eval]
ax.plot(t_eval, analytical_values, 'k--', linewidth=3, label='Analytical Solution', alpha=0.8)
# Plot numerical solution
colors = ['red', 'blue', 'green', 'orange', 'purple']
for i, (method_name, method_data) in enumerate(results[ode_name].items()):
if i >= len(colors):
break
sol = method_data['solution']
if ode_name == 'harmonic':
y_values = sol.y[0, :] # Position component
else:
y_values = sol.y[0, :]
ax.plot(sol.t, y_values, color=colors[i], linewidth=2,
label=f'{method_name} (Error: {method_data["error"]:.2e})', alpha=0.7)
ax.set_title(f'{description}')
ax.set_xlabel('Time t')
ax.set_ylabel('y(t)')
ax.legend()
ax.grid(True, alpha=0.3)
plot_idx += 1
# Phase space plot (only for harmonic oscillator)
if 'harmonic' in results:
ax = axes[3]
for method_name, method_data in list(results['harmonic'].items())[:3]: # Only show first 3 methods
sol = method_data['solution']
y_pos = sol.y[0, :]
y_vel = sol.y[1, :]
ax.plot(y_pos, y_vel, linewidth=2, label=f'{method_name}', alpha=0.7)
# Phase space trajectory of analytical solution
y_pos_analytical = [harmonic_analytical(t) for t in t_eval]
y_vel_analytical = [-2 * harmonic_analytical(t, y0=1, v0=0, omega=2) *
np.sin(2*t) for t in t_eval] # Simplified velocity
ax.plot(y_pos_analytical, y_vel_analytical, 'k--', linewidth=3,
label='Analytical Solution', alpha=0.8)
ax.set_title('Harmonic Oscillator Phase Space Trajectory')
ax.set_xlabel('Position y')
ax.set_ylabel('Velocity dy/dt')
ax.legend()
ax.grid(True, alpha=0.3)
ax.axis('equal')
plt.tight_layout()
plt.show()2. Stiff Equations and Higher-Order Equations
python
# Stiff equation examples
def stiff_ode_example():
"""Stiff equation examples"""
# Van der Pol oscillator (stiff parameter μ)
def van_der_pol(t, y, mu=10):
y1, y2 = y
return [y2, mu * (1 - y1**2) * y2 - y1]
# Chemical reaction kinetics (Robertson problem)
def robertson(t, y):
y1, y2, y3 = y
return [-0.04*y1 + 1e4*y2*y3,
0.04*y1 - 1e4*y2*y3 - 3e7*y2**2,
3e7*y2**2]
return {
'van_der_pol': (van_der_pol, [2, 0], 'Van der Pol Oscillator'),
'robertson': (robertson, [1, 0, 0], 'Robertson Chemical Reaction')
}
# Get stiff equation examples
stiff_examples = stiff_ode_example()
# Compare methods suitable for stiff equations
stiff_methods = {
'Radau': 'Radau', # Implicit method, suitable for stiff equations
'BDF': 'BDF', # Backward differentiation, suitable for stiff equations
'RK45': 'RK45' # Explicit method, for comparison
}
print("Stiff Equation Solving Comparison:")
print(f"{'Equation':>15} {'Method':>10} {'Success':>8} {'Steps':>8} {'Time(ms)':>10} {'Final Time':>10}")
print("-" * 70)
t_span_stiff = (0, 20)
t_eval_stiff = np.linspace(0, 20, 200)
stiff_results = {}
for eq_name, (eq_func, y0, description) in stiff_examples.items():
stiff_results[eq_name] = {}
for method_name, method in stiff_methods.items():
try:
start_time = time.time()
if eq_name == 'van_der_pol':
sol = solve_ivp(lambda t, y: van_der_pol(t, y, mu=10),
t_span_stiff, y0, method=method,
dense_output=True, rtol=1e-6)
else:
sol = solve_ivp(eq_func, t_span_stiff, y0, method=method,
dense_output=True, rtol=1e-6)
exec_time = (time.time() - start_time) * 1000
success = "Yes" if sol.success else "No"
n_steps = len(sol.t)
final_time = sol.t[-1] if len(sol.t) > 0 else 0
stiff_results[eq_name][method_name] = {
'solution': sol,
'success': sol.success,
'n_steps': n_steps,
'time': exec_time
}
print(f"{eq_name:>15} {method_name:>10} {success:>8} {n_steps:>8} {exec_time:>10.2f} {final_time:>10.2f}")
except Exception as e:
print(f"{eq_name:>15} {method_name:>10} {'No':>8} {'N/A':>8} {'N/A':>10} {'N/A':>10}")
stiff_results[eq_name][method_name] = {'success': False}
print()
# Visualize stiff equation solutions
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
# Van der Pol oscillator
if 'van_der_pol' in stiff_results:
# Time series
ax1 = axes[0, 0]
ax2 = axes[0, 1]
for method_name, result in stiff_results['van_der_pol'].items():
if result['success']:
sol = result['solution']
t_dense = np.linspace(0, 20, 1000)
y_dense = sol.sol(t_dense)
ax1.plot(t_dense, y_dense[0], label=f'{method_name}', linewidth=2)
ax1.set_title('Van der Pol Oscillator - Time Series')
ax1.set_xlabel('Time t')
ax1.set_ylabel('y1(t)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Phase space
for method_name, result in stiff_results['van_der_pol'].items():
if result['success']:
sol = result['solution']
t_dense = np.linspace(0, 20, 1000)
y_dense = sol.sol(t_dense)
ax2.plot(y_dense[0], y_dense[1], label=f'{method_name}', linewidth=2)
ax2.set_title('Van der Pol Oscillator - Phase Space')
ax2.set_xlabel('y1')
ax2.set_ylabel('y2')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Robertson reaction
if 'robertson' in stiff_results:
ax3 = axes[1, 0]
ax4 = axes[1, 1]
# Concentration over time
for method_name, result in stiff_results['robertson'].items():
if result['success']:
sol = result['solution']
t_dense = np.linspace(0, 20, 1000)
y_dense = sol.sol(t_dense)
ax3.semilogy(t_dense, y_dense[0], label=f'{method_name} - A', linewidth=2)
ax3.semilogy(t_dense, y_dense[1], '--', label=f'{method_name} - B', linewidth=2)
ax3.semilogy(t_dense, y_dense[2], ':', label=f'{method_name} - C', linewidth=2)
ax3.set_title('Robertson Reaction - Concentration Changes')
ax3.set_xlabel('Time t')
ax3.set_ylabel('Concentration (Log Scale)')
ax3.legend()
ax3.grid(True, alpha=0.3)
# Conservation quantity check
for method_name, result in stiff_results['robertson'].items():
if result['success']:
sol = result['solution']
t_dense = np.linspace(0, 20, 1000)
y_dense = sol.sol(t_dense)
conservation = y_dense[0] + y_dense[1] + y_dense[2]
ax4.plot(t_dense, conservation - 1, label=f'{method_name}', linewidth=2)
ax4.set_title('Robertson Reaction - Conservation Deviation')
ax4.set_xlabel('Time t')
ax4.set_ylabel('Total Concentration - 1')
ax4.legend()
ax4.grid(True, alpha=0.3)
ax4.ticklabel_format(style='scientific', axis='y', scilimits=(0,0))
plt.tight_layout()
plt.show()3. Boundary Value Problems
python
# Boundary value problem solving
def boundary_value_problems():
"""Boundary value problem examples"""
# Example 1: Second-order linear BVP
# y'' + y = 0, y(0) = 0, y(π) = 0
# Analytical solution: y(x) = A*sin(x), from boundary conditions A = 0 or sin(π) = 0
def linear_bvp(x, y):
return np.vstack((y[1], -y[0]))
def linear_bc(ya, yb):
return np.array([ya[0], yb[0]]) # y(0) = 0, y(π) = 0
# Example 2: Nonlinear BVP (Bratu problem)
# y'' + λ*exp(y) = 0, y(0) = y(1) = 0
def bratu_bvp(x, y, lam=1):
return np.vstack((y[1], -lam * np.exp(y[0])))
def bratu_bc(ya, yb):
return np.array([ya[0], yb[0]]) # y(0) = y(1) = 0
return {
'linear': (linear_bvp, linear_bc, [0, np.pi], 'Linear BVP: y\'\' + y = 0'),
'bratu': (bratu_bvp, bratu_bc, [0, 1], 'Bratu Problem: y\'\' + λe^y = 0')
}
# Get BVP examples
bvp_examples = boundary_value_problems()
print("Boundary Value Problem Solving:")
print(f"{'Problem':>15} {'Convergence':>8} {'Residual':>12} {'Iterations':>10}")
print("-" * 50)
bvp_results = {}
for bvp_name, (bvp_func, bc_func, domain, description) in bvp_examples.items():
try:
# Initial grid
x = np.linspace(domain[0], domain[1], 11)
# Initial guess
if bvp_name == 'linear':
y_guess = np.zeros((2, x.size))
y_guess[0] = np.sin(x) # Guess solution shape
else: # bratu
y_guess = np.zeros((2, x.size))
y_guess[0] = x * (1 - x) # Initial guess satisfying boundary conditions
# Solve BVP
if bvp_name == 'bratu':
sol = solve_bvp(lambda x, y: bratu_bvp(x, y, lam=1), bc_func, x, y_guess)
else:
sol = solve_bvp(bvp_func, bc_func, x, y_guess)
converged = "Yes" if sol.success else "No"
residual = np.max(np.abs(sol.rms_residuals)) if sol.success else float('inf')
n_iter = sol.niter if hasattr(sol, 'niter') else 'N/A'
bvp_results[bvp_name] = {
'solution': sol,
'success': sol.success,
'residual': residual
}
print(f"{bvp_name:>15} {converged:>8} {residual:>12.2e} {n_iter:>10}")
except Exception as e:
print(f"{bvp_name:>15} {'No':>8} {'N/A':>12} {'N/A':>10}")
bvp_results[bvp_name] = {'success': False}
# Visualize BVP solutions
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
for i, (bvp_name, result) in enumerate(bvp_results.items()):
if i >= 2 or not result['success']:
continue
ax = axes[i]
sol = result['solution']
# Plot solution
x_plot = np.linspace(sol.x[0], sol.x[-1], 200)
y_plot = sol.sol(x_plot)
ax.plot(x_plot, y_plot[0], 'b-', linewidth=3, label='Numerical Solution')
ax.plot(sol.x, sol.y[0], 'ro', markersize=6, label='Grid Points')
# Add boundary condition markers
ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
ax.plot([sol.x[0], sol.x[-1]], [0, 0], 'go', markersize=8, label='Boundary Conditions')
ax.set_title(f'{bvp_examples[bvp_name][3]}')
ax.set_xlabel('x')
ax.set_ylabel('y(x)')
ax.legend()
ax.grid(True, alpha=0.3)
# Add residual information
ax.text(0.05, 0.95, f'Max Residual: {result["residual"]:.2e}',
transform=ax.transAxes, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
plt.tight_layout()
plt.show()Practical Application Cases
1. Population Dynamics Modeling
python
# Population dynamics modeling
def population_dynamics():
"""Population dynamics models"""
# Lotka-Volterra predator-prey model
def lotka_volterra(t, y, a=1, b=0.1, c=1.5, d=0.075):
"""Lotka-Volterra equations
y[0]: Prey population
y[1]: Predator population
"""
prey, predator = y
dprey_dt = a * prey - b * prey * predator
dpredator_dt = -c * predator + d * prey * predator
return [dprey_dt, dpredator_dt]
# Competition model
def competition_model(t, y, r1=1, r2=0.8, K1=100, K2=80, a12=0.5, a21=0.7):
"""Interspecific competition model"""
N1, N2 = y
dN1_dt = r1 * N1 * (1 - (N1 + a12 * N2) / K1)
dN2_dt = r2 * N2 * (1 - (N2 + a21 * N1) / K2)
return [dN1_dt, dN2_dt]
# SIR epidemic model
def sir_model(t, y, beta=0.3, gamma=0.1):
"""SIR epidemic model"""
S, I, R = y
N = S + I + R
dS_dt = -beta * S * I / N
dI_dt = beta * S * I / N - gamma * I
dR_dt = gamma * I
return [dS_dt, dI_dt, dR_dt]
return {
'lotka_volterra': (lotka_volterra, [10, 5], 'Lotka-Volterra Predator-Prey Model'),
'competition': (competition_model, [50, 40], 'Interspecific Competition Model'),
'sir': (sir_model, [990, 10, 0], 'SIR Epidemic Model')
}
# Get population dynamics models
population_models = population_dynamics()
# Solve population dynamics equations
t_span = (0, 20)
t_eval = np.linspace(0, 20, 1000)
print("Population Dynamics Model Solving:")
population_results = {}
for model_name, (model_func, y0, description) in population_models.items():
try:
sol = solve_ivp(model_func, t_span, y0, t_eval=t_eval, method='RK45')
population_results[model_name] = {
'solution': sol,
'success': sol.success,
'description': description
}
print(f"{model_name:>15}: {'Success' if sol.success else 'Failed'}")
except Exception as e:
print(f"{model_name:>15}: Error - {e}")
population_results[model_name] = {'success': False}
# Visualize population dynamics results
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes = axes.flatten()
plot_idx = 0
for model_name, result in population_results.items():
if not result['success'] or plot_idx >= 6:
continue
sol = result['solution']
description = result['description']
# Time series plot
ax1 = axes[plot_idx]
plot_idx += 1
if model_name == 'lotka_volterra':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Prey')
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Predator')
ax1.set_ylabel('Population')
elif model_name == 'competition':
ax1.plot(sol.t, sol.y[0], 'g-', linewidth=2, label='Species 1')
ax1.plot(sol.t, sol.y[1], 'm-', linewidth=2, label='Species 2')
ax1.set_ylabel('Population')
elif model_name == 'sir':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Susceptible (S)')
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Infected (I)')
ax1.plot(sol.t, sol.y[2], 'g-', linewidth=2, label='Recovered (R)')
ax1.set_ylabel('Population')
ax1.set_title(f'{description} - Time Series')
ax1.set_xlabel('Time')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Phase space plot (except SIR model)
if model_name != 'sir' and plot_idx < 6:
ax2 = axes[plot_idx]
plot_idx += 1
ax2.plot(sol.y[0], sol.y[1], 'purple', linewidth=2)
ax2.plot(sol.y[0, 0], sol.y[1, 0], 'go', markersize=8, label='Start')
ax2.plot(sol.y[0, -1], sol.y[1, -1], 'ro', markersize=8, label='End')
if model_name == 'lotka_volterra':
ax2.set_xlabel('Prey Population')
ax2.set_ylabel('Predator Population')
ax2.set_title('Lotka-Volterra Phase Space Trajectory')
elif model_name == 'competition':
ax2.set_xlabel('Species 1 Population')
ax2.set_ylabel('Species 2 Population')
ax2.set_title('Competition Model Phase Space Trajectory')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Special analysis for SIR model
if 'sir' in population_results and population_results['sir']['success']:
if plot_idx < 6:
ax_sir = axes[plot_idx]
sol_sir = population_results['sir']['solution']
# Calculate basic reproduction number R0 and peak time
I_max_idx = np.argmax(sol_sir.y[1])
t_peak = sol_sir.t[I_max_idx]
I_peak = sol_sir.y[1, I_max_idx]
# Plot infected curve and peak
ax_sir.plot(sol_sir.t, sol_sir.y[1], 'r-', linewidth=3, label='Infected')
ax_sir.axvline(x=t_peak, color='orange', linestyle='--',
label=f'Peak Time: {t_peak:.1f}')
ax_sir.axhline(y=I_peak, color='orange', linestyle='--',
label=f'Peak Population: {I_peak:.0f}')
ax_sir.set_title('SIR Model - Infected Peak Analysis')
ax_sir.set_xlabel('Time')
ax_sir.set_ylabel('Infected Population')
ax_sir.legend()
ax_sir.grid(True, alpha=0.3)
# Add statistical information
total_infected = 1000 - sol_sir.y[0, -1] # Final susceptible reduction
ax_sir.text(0.05, 0.95, f'Total Infected: {total_infected:.0f}\nInfection Rate: {total_infected/1000*100:.1f}%',
transform=ax_sir.transAxes, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
plt.tight_layout()
plt.show()
# Parameter sensitivity analysis
print("\nParameter Sensitivity Analysis (SIR Model):")
if 'sir' in population_results and population_results['sir']['success']:
beta_values = [0.1, 0.2, 0.3, 0.4, 0.5]
gamma_values = [0.05, 0.1, 0.15, 0.2]
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
# β parameter effect
for beta in beta_values:
sol = solve_ivp(lambda t, y: sir_model(t, y, beta=beta, gamma=0.1),
t_span, [990, 10, 0], t_eval=t_eval)
if sol.success:
axes[0].plot(sol.t, sol.y[1], linewidth=2, label=f'β={beta}')
axes[0].set_title('Effect of Transmission Rate β on Infected Peak')
axes[0].set_xlabel('Time')
axes[0].set_ylabel('Infected Population')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# γ parameter effect
for gamma in gamma_values:
sol = solve_ivp(lambda t, y: sir_model(t, y, beta=0.3, gamma=gamma),
t_span, [990, 10, 0], t_eval=t_eval)
if sol.success:
axes[1].plot(sol.t, sol.y[1], linewidth=2, label=f'γ={gamma}')
axes[1].set_title('Effect of Recovery Rate γ on Infected Peak')
axes[1].set_xlabel('Time')
axes[1].set_ylabel('Infected Population')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()2. Physical System Modeling
python
# Physical system modeling
def physics_systems():
"""Physical system differential equations"""
# Pendulum equation
def pendulum(t, y, g=9.81, L=1.0, damping=0.1):
"""Damped pendulum equation
θ'' + (damping/m)θ' + (g/L)sin(θ) = 0
"""
theta, omega = y
return [omega, -damping * omega - (g/L) * np.sin(theta)]
# Spring-mass-damper system
def spring_mass_damper(t, y, m=1.0, k=4.0, c=0.5):
"""Spring-mass-damper system
mx'' + cx' + kx = 0
"""
x, v = y
return [v, -(c/m) * v - (k/m) * x]
# RLC circuit
def rlc_circuit(t, y, R=1.0, L=1.0, C=1.0, V0=0):
"""RLC circuit equation
L(dI/dt) + R*I + Q/C = V0
dQ/dt = I
"""
Q, I = y # Charge and current
return [I, (V0 - R*I - Q/C) / L]
# Lorenz equations (chaotic system)
def lorenz(t, y, sigma=10, rho=28, beta=8/3):
"""Lorenz equations"""
x, y_coord, z = y
return [sigma * (y_coord - x),
x * (rho - z) - y_coord,
x * y_coord - beta * z]
return {
'pendulum': (pendulum, [np.pi/4, 0], 'Damped Pendulum'),
'spring_mass': (spring_mass_damper, [1, 0], 'Spring-Mass-Damper System'),
'rlc_circuit': (rlc_circuit, [0, 0], 'RLC Circuit'),
'lorenz': (lorenz, [1, 1, 1], 'Lorenz Chaotic System')
}
# Get physical system models
physics_models = physics_systems()
# Solve physical system equations
t_span_physics = (0, 10)
t_eval_physics = np.linspace(0, 10, 2000)
print("Physical System Modeling:")
physics_results = {}
for model_name, (model_func, y0, description) in physics_models.items():
try:
if model_name == 'lorenz':
# Lorenz system needs longer time
t_span_lorenz = (0, 30)
t_eval_lorenz = np.linspace(0, 30, 5000)
sol = solve_ivp(model_func, t_span_lorenz, y0, t_eval=t_eval_lorenz, method='RK45')
else:
sol = solve_ivp(model_func, t_span_physics, y0, t_eval=t_eval_physics, method='RK45')
physics_results[model_name] = {
'solution': sol,
'success': sol.success,
'description': description
}
print(f"{model_name:>15}: {'Success' if sol.success else 'Failed'}")
except Exception as e:
print(f"{model_name:>15}: Error - {e}")
physics_results[model_name] = {'success': False}
# Visualize physical system results
fig, axes = plt.subplots(2, 4, figsize=(20, 10))
axes = axes.flatten()
plot_idx = 0
for model_name, result in physics_results.items():
if not result['success']:
continue
sol = result['solution']
description = result['description']
# Time series plot
ax1 = axes[plot_idx]
plot_idx += 1
if model_name == 'pendulum':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Angle θ')
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Angular Velocity ω')
ax1.set_ylabel('Angle/Angular Velocity')
elif model_name == 'spring_mass':
ax1.plot(sol.t, sol.y[0], 'g-', linewidth=2, label='Displacement x')
ax1.plot(sol.t, sol.y[1], 'm-', linewidth=2, label='Velocity v')
ax1.set_ylabel('Displacement/Velocity')
elif model_name == 'rlc_circuit':
ax1.plot(sol.t, sol.y[0], 'c-', linewidth=2, label='Charge Q')
ax1.plot(sol.t, sol.y[1], 'orange', linewidth=2, label='Current I')
ax1.set_ylabel('Charge/Current')
elif model_name == 'lorenz':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=1, label='x', alpha=0.8)
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=1, label='y', alpha=0.8)
ax1.plot(sol.t, sol.y[2], 'g-', linewidth=1, label='z', alpha=0.8)
ax1.set_ylabel('Coordinate Value')
ax1.set_title(f'{description} - Time Series')
ax1.set_xlabel('Time')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Phase space plot
if plot_idx < 8:
ax2 = axes[plot_idx]
plot_idx += 1
if model_name == 'lorenz':
# 3D trajectory projected to 2D
ax2.plot(sol.y[0], sol.y[1], 'purple', linewidth=0.5, alpha=0.7)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('Lorenz Attractor (x-y Projection)')
else:
ax2.plot(sol.y[0], sol.y[1], 'purple', linewidth=2)
ax2.plot(sol.y[0, 0], sol.y[1, 0], 'go', markersize=8, label='Start')
if model_name == 'pendulum':
ax2.set_xlabel('Angle θ')
ax2.set_ylabel('Angular Velocity ω')
ax2.set_title('Pendulum Phase Space Trajectory')
elif model_name == 'spring_mass':
ax2.set_xlabel('Displacement x')
ax2.set_ylabel('Velocity v')
ax2.set_title('Spring System Phase Space Trajectory')
elif model_name == 'rlc_circuit':
ax2.set_xlabel('Charge Q')
ax2.set_ylabel('Current I')
ax2.set_title('RLC Circuit Phase Space Trajectory')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 3D visualization of Lorenz system
if 'lorenz' in physics_results and physics_results['lorenz']['success']:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(15, 5))
sol_lorenz = physics_results['lorenz']['solution']
# 3D trajectory
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot(sol_lorenz.y[0], sol_lorenz.y[1], sol_lorenz.y[2],
'b-', linewidth=0.5, alpha=0.8)
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')
ax1.set_title('Lorenz Attractor 3D Trajectory')
# x-z projection
ax2 = fig.add_subplot(132)
ax2.plot(sol_lorenz.y[0], sol_lorenz.y[2], 'r-', linewidth=0.5, alpha=0.8)
ax2.set_xlabel('X')
ax2.set_ylabel('Z')
ax2.set_title('Lorenz Attractor (x-z Projection)')
ax2.grid(True, alpha=0.3)
# y-z projection
ax3 = fig.add_subplot(133)
ax3.plot(sol_lorenz.y[1], sol_lorenz.y[2], 'g-', linewidth=0.5, alpha=0.8)
ax3.set_xlabel('Y')
ax3.set_ylabel('Z')
ax3.set_title('Lorenz Attractor (y-z Projection)')
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
## Performance Optimization Tips
### 1. Algorithm Selection and Parameter Tuning
```python
# Performance optimization examples
def performance_optimization():
"""Integration and ODE solving performance optimization"""
# Test functions
def test_integration_performance():
"""Integration performance test"""
def smooth_func(x):
return np.exp(-x**2) * np.sin(x)
def oscillatory_func(x):
return np.sin(100*x) / (1 + x**2)
methods = {
'quad': lambda f, a, b: quad(f, a, b),
'fixed_quad_5': lambda f, a, b: fixed_quad(f, a, b, n=5),
'fixed_quad_10': lambda f, a, b: fixed_quad(f, a, b, n=10),
'simpson_1000': lambda f, a, b: (simpson([f(x) for x in np.linspace(a, b, 1001)], dx=(b-a)/1000), 0)
}
functions = {'smooth': smooth_func, 'oscillatory': oscillatory_func}
print("Integration Method Performance Comparison:")
print(f"{'Function':>12} {'Method':>15} {'Result':>12} {'Time(ms)':>10} {'Precision':>12}")
print("-" * 70)
for func_name, func in functions.items():
# Reference result
ref_result, _ = quad(func, 0, 1)
for method_name, method in methods.items():
times = []
for _ in range(10): # Multiple tests for averaging
start = time.time()
result, _ = method(func, 0, 1)
times.append((time.time() - start) * 1000)
avg_time = np.mean(times)
accuracy = abs(result - ref_result)
print(f"{func_name:>12} {method_name:>15} {result:>12.6f} {avg_time:>10.2f} {accuracy:>12.2e}")
print()
def test_ode_performance():
"""ODE solving performance test"""
def stiff_ode(t, y):
return [-1000*y[0] + y[1], y[0] - y[1]]
def nonstiff_ode(t, y):
return [-y[0] + y[1], y[0] - y[1]]
methods = ['RK45', 'RK23', 'Radau', 'BDF']
odes = {
'non_stiff': (nonstiff_ode, [1, 0]),
'stiff': (stiff_ode, [1, 0])
}
print("ODE Solving Method Performance Comparison:")
print(f"{'Equation Type':>12} {'Method':>10} {'Success':>6} {'Steps':>8} {'Time(ms)':>10} {'Precision':>12}")
print("-" * 70)
for ode_name, (ode_func, y0) in odes.items():
for method in methods:
try:
start = time.time()
sol = solve_ivp(ode_func, (0, 1), y0, method=method, rtol=1e-6)
exec_time = (time.time() - start) * 1000
success = "Yes" if sol.success else "No"
n_steps = len(sol.t)
# Simple precision estimation (compared with most accurate method)
if method == 'Radau':
reference_sol = sol
accuracy = 0
else:
try:
ref_sol = solve_ivp(ode_func, (0, 1), y0, method='Radau', rtol=1e-8)
if ref_sol.success:
accuracy = np.max(np.abs(sol.y[:, -1] - ref_sol.y[:, -1]))
else:
accuracy = float('inf')
except:
accuracy = float('inf')
print(f"{ode_name:>12} {method:>10} {success:>6} {n_steps:>8} {exec_time:>10.2f} {accuracy:>12.2e}")
except Exception as e:
print(f"{ode_name:>12} {method:>10} {'No':>6} {'N/A':>8} {'N/A':>10} {'N/A':>12}")
print()
# Execute performance tests
test_integration_performance()
test_ode_performance()
# Execute performance optimization tests
performance_optimization()2. Memory Optimization and Parallel Computing
python
# Memory optimization and parallel computing
def memory_and_parallel_optimization():
"""Memory optimization and parallel computing examples"""
# Memory optimization for large-scale integration
def memory_efficient_integration():
"""Memory-efficient integration calculation"""
def large_array_function(x):
# Function requiring large memory
return np.sum(np.sin(np.arange(1000) * x))
# Chunked integration
def chunked_integration(func, a, b, n_chunks=10):
"""Chunked integration to save memory"""
chunk_size = (b - a) / n_chunks
total_integral = 0
for i in range(n_chunks):
chunk_a = a + i * chunk_size
chunk_b = a + (i + 1) * chunk_size
chunk_integral, _ = quad(func, chunk_a, chunk_b)
total_integral += chunk_integral
return total_integral
# Compare memory usage
print("Memory Optimization Integration Comparison:")
# Direct integration
start_time = time.time()
direct_result, _ = quad(large_array_function, 0, 1)
direct_time = time.time() - start_time
# Chunked integration
start_time = time.time()
chunked_result = chunked_integration(large_array_function, 0, 1, n_chunks=5)
chunked_time = time.time() - start_time
print(f"Direct Integration: {direct_result:.6f}, Time: {direct_time*1000:.2f}ms")
print(f"Chunked Integration: {chunked_result:.6f}, Time: {chunked_time*1000:.2f}ms")
print(f"Relative Error: {abs(direct_result - chunked_result)/abs(direct_result):.2e}")
print()
# Parallel ODE solving
def parallel_ode_solving():
"""Parallel ODE solving examples"""
def parameter_study_ode(t, y, param):
"""Parameterized ODE"""
return [-param * y[0], y[0] - y[1]]
# Parameter range
parameters = np.linspace(0.1, 10, 20)
# Serial solving
start_time = time.time()
serial_results = []
for param in parameters:
sol = solve_ivp(lambda t, y: parameter_study_ode(t, y, param),
(0, 5), [1, 0], method='RK45')
if sol.success:
serial_results.append(sol.y[0, -1])
else:
serial_results.append(np.nan)
serial_time = time.time() - start_time
print(f"Serial Solving Time: {serial_time*1000:.2f}ms")
print(f"Successful Solving: {np.sum(~np.isnan(serial_results))}/{len(parameters)}")
return serial_results, parameters
# Execute optimization tests
memory_efficient_integration()
results, params = parallel_ode_solving()
# Visualize parameter study results
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(params, results, 'bo-', linewidth=2, markersize=6)
plt.xlabel('Parameter Value')
plt.ylabel('Final State y[0](T)')
plt.title('Parameter Sensitivity Analysis')
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
# Show convergence
valid_results = np.array(results)[~np.isnan(results)]
valid_params = params[~np.isnan(results)]
if len(valid_results) > 1:
# Calculate gradient
gradient = np.gradient(valid_results, valid_params)
plt.plot(valid_params, gradient, 'ro-', linewidth=2, markersize=6)
plt.xlabel('Parameter Value')
plt.ylabel('Sensitivity dy/dp')
plt.title('Parameter Sensitivity Gradient')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Execute memory and parallel optimization
memory_and_parallel_optimization()
## Summary
SciPy's integration and differential equations module provides powerful and flexible numerical computing tools:
### Integration Features
- **One-Dimensional Integration**: `quad`, `simpson`, `trapz` and other methods suitable for different precision requirements
- **Multiple Integration**: `dblquad`, `tplquad`, `nquad` support complex region integration
- **Special Integration**: Handle singularities, infinite intervals, oscillatory functions, and other special cases
### Differential Equation Solving
- **Initial Value Problems**: `solve_ivp` provides multiple high-precision solvers
- **Boundary Value Problems**: `solve_bvp` handles complex boundary conditions
- **Stiff Equations**: `Radau`, `BDF` and other implicit methods handle stiff systems
### Application Fields
- **Scientific Computing**: Physical modeling, engineering analysis, numerical simulation
- **Mathematical Biology**: Population dynamics, epidemic models, ecosystems
- **Engineering Applications**: Control systems, signal processing, circuit analysis
### Performance Optimization
- **Algorithm Selection**: Choose appropriate solvers based on problem characteristics
- **Parameter Tuning**: Reasonably set tolerances, step sizes, and other parameters
- **Memory Management**: Chunked processing for large-scale problems