积分和微分方程
SciPy 提供了强大的数值积分和微分方程求解功能，主要通过 scipy.integrate 模块实现。该模块包含了各种数值积分方法、常微分方程(ODE)求解器、偏微分方程(PDE)求解工具等。这些功能在科学计算、工程分析、物理建模等领域有着广泛的应用。本章将详细介绍如何使用 SciPy 进行数值积分和微分方程求解。
scipy.integrate 模块概述
scipy.integrate 模块包含以下主要功能：
- 数值积分（定积分和不定积分）
- 常微分方程求解
- 边值问题求解
- 积分方程求解
- 多重积分计算
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')
# 设置绘图样式
plt.style.use('seaborn-v0_8')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 查看 integrate 模块的主要功能
print("scipy.integrate 主要功能:")
functions = [attr for attr in dir(integrate) if not attr.startswith('_')]
print(f"总共 {len(functions)} 个函数和类")
print("积分函数:", [f for f in functions if 'quad' in f or 'int' in f][:10])
print("ODE求解器:", [f for f in functions if 'ode' in f or 'solve' in f][:10])数值积分
1. 一维积分
python
# 定义测试函数
def test_functions():
"""定义各种测试函数"""
functions = {
'多项式': lambda x: x**3 - 2*x**2 + x - 1,
'三角函数': lambda x: np.sin(x) * np.cos(x),
'指数函数': lambda x: np.exp(-x**2),
'有理函数': lambda x: 1 / (1 + x**2),
'振荡函数': lambda x: np.sin(10*x) / x if x != 0 else 10,
'不连续函数': lambda x: np.where(x < 0.5, x**2, 2*x - 1)
}
# 解析解（用于验证）
analytical_solutions = {
'多项式': 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),
'三角函数': lambda a, b: (np.sin(b)**2 - np.sin(a)**2) / 2,
'有理函数': lambda a, b: np.arctan(b) - np.arctan(a)
}
return functions, analytical_solutions
# 获取测试函数
test_funcs, analytical_sols = test_functions()
# 使用不同的积分方法
integration_methods = {
'quad (自适应)': lambda f, a, b: quad(f, a, b),
'fixed_quad (高斯)': lambda f, a, b: fixed_quad(f, a, b, n=5),
'quadrature (自适应高斯)': lambda f, a, b: quadrature(f, a, b)
}
# 积分区间
a, b = 0, 1
print(f"数值积分比较 (区间 [{a}, {b}]):")
print(f"{'函数':>15} {'方法':>15} {'数值结果':>12} {'误差估计':>12} {'解析解':>12} {'绝对误差':>12}")
print("-" * 90)
for func_name, func in test_funcs.items():
# 计算解析解（如果有）
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)
# 计算与解析解的误差
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} {'失败':>12} {'N/A':>12} {'N/A':>12} {'N/A':>12}")
print()
# 可视化积分过程
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
# 计算函数值
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()
# 添加积分值标注
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'绘图失败\n{str(e)}', ha='center', va='center')
axes[i].set_title(f'{name} (绘图失败)')
plt.tight_layout()
plt.show()2. 多重积分
python
# 二重积分示例
def double_integration_examples():
"""二重积分示例"""
# 定义二元函数
functions_2d = {
'简单多项式': lambda y, x: x*y,
'高斯函数': lambda y, x: np.exp(-(x**2 + y**2)),
'三角函数': lambda y, x: np.sin(x) * np.cos(y),
'有理函数': lambda y, x: 1 / (1 + x**2 + y**2)
}
# 积分区域
regions = {
'矩形区域': {'x_bounds': [0, 1], 'y_bounds': [0, 1]},
'三角形区域': {'x_bounds': [0, 1], 'y_bounds': lambda x: [0, x]},
'圆形区域': {'x_bounds': [-1, 1], 'y_bounds': lambda x: [-np.sqrt(1-x**2), np.sqrt(1-x**2)]}
}
print("二重积分计算:")
print(f"{'函数':>15} {'区域':>15} {'积分值':>12} {'误差估计':>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']):
# 变限积分
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:
# 常限积分
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} {'失败':>12} {'N/A':>12}")
results[func_name][region_name] = (None, None)
return results
# 执行二重积分计算
double_integral_results = double_integration_examples()
# 三重积分示例
print("\n三重积分计算:")
# 定义三元函数
def triple_function(z, y, x):
return x*y*z
# 计算三重积分
try:
# 在单位立方体上积分
result_3d, error_3d = tplquad(triple_function, 0, 1, 0, 1, 0, 1)
print(f"∭xyz dxdydz (单位立方体) = {result_3d:.6f} ± {error_3d:.2e}")
# 解析解: 1/8
analytical_3d = 1/8
print(f"解析解 = {analytical_3d:.6f}")
print(f"绝对误差 = {abs(result_3d - analytical_3d):.2e}")
except Exception as e:
print(f"三重积分计算失败: {e}")
# 可视化二重积分
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()
# 创建网格
x = np.linspace(-1, 1, 50)
y = np.linspace(-1, 1, 50)
X, Y = np.meshgrid(x, y)
# 可视化不同的二元函数
functions_for_plot = {
'高斯函数': lambda x, y: np.exp(-(x**2 + y**2)),
'三角函数': lambda x, y: np.sin(x) * np.cos(y),
'有理函数': lambda x, y: 1 / (1 + x**2 + y**2),
'波函数': 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])
# 添加积分区域标记
circle = plt.Circle((0, 0), 1, fill=False, color='red', linewidth=2, linestyle='--')
axes[i].add_patch(circle)
axes[i].text(0, -1.3, '红色虚线: 圆形积分区域', ha='center', color='red')
plt.tight_layout()
plt.show()3. 数值积分方法比较
python
# 比较不同数值积分方法的精度和效率
import time
def compare_integration_methods():
"""比较不同积分方法"""
# 测试函数
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,
'振荡函数': oscillatory_function,
'奇异函数': singular_function
}
# 积分方法
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)
}
# 积分区间
intervals = {
'光滑函数': (0, 2),
'振荡函数': (0, 1),
'奇异函数': (0.001, 1) # 避免奇点
}
print("积分方法比较:")
print(f"{'函数':>12} {'方法':>12} {'结果':>12} {'时间(ms)':>10} {'相对误差':>12}")
print("-" * 70)
for func_name, func in test_functions.items():
a, b = intervals[func_name]
# 使用quad作为参考
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} {'失败':>12} {'N/A':>10} {'N/A':>12}")
print()
# 执行方法比较
compare_integration_methods()
# 可视化积分精度
def visualize_integration_accuracy():
"""可视化积分精度"""
def test_func(x):
return np.sin(x) * np.exp(-x/2)
# 解析解
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)
# 不同的网格密度
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法
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)
# 梯形法
trapz_result = trapz(y, dx=dx)
methods_accuracy['Trapezoid'].append(abs(trapz_result - true_value))
# Quad法（参考）
quad_result, _ = quad(test_func, a, b)
methods_accuracy['Quad'].append(abs(quad_result - true_value))
# 绘制精度比较
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('被积函数')
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误差太小，单独显示
plt.loglog(n_points, errors, 'o-', label=method, linewidth=2, markersize=6)
plt.xlabel('网格点数')
plt.ylabel('绝对误差')
plt.title('积分精度比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(2, 2, 3)
# 显示收敛阶
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收敛阶: {-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收敛阶: {-trapz_order:.2f}', transform=plt.gca().transAxes)
plt.text(0.1, 0.6, f'真实积分值: {true_value:.6f}', transform=plt.gca().transAxes)
plt.text(0.1, 0.5, f'Quad结果: {quad_result:.6f}', transform=plt.gca().transAxes)
plt.text(0.1, 0.4, f'Quad误差: {abs(quad_result - true_value):.2e}', transform=plt.gca().transAxes)
plt.title('收敛性分析')
plt.axis('off')
plt.subplot(2, 2, 4)
# 显示不同方法的最终精度
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('绝对误差')
plt.title('最终精度比较')
plt.grid(True, alpha=0.3)
# 添加数值标签
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()
# 执行精度可视化
visualize_integration_accuracy()常微分方程求解
1. 初值问题
python
# 常微分方程初值问题求解
def ode_examples():
"""常微分方程示例"""
# 示例1: 一阶线性ODE
# dy/dt = -2y + 1, y(0) = 0
# 解析解: 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))
# 示例2: 非线性ODE (Logistic方程)
# 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))
# 示例3: 振荡器 (简谐振动)
# 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')
}
# 获取ODE示例
ode_examples = ode_examples()
# 求解ODE并比较方法
t_span = (0, 5)
t_eval = np.linspace(0, 5, 100)
print("常微分方程求解比较:")
print(f"{'方程':>15} {'方法':>15} {'最终值':>12} {'解析解':>12} {'绝对误差':>12} {'时间(ms)':>10}")
print("-" * 85)
# 不同的求解方法
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] = {}
# 计算解析解
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
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] # 位置分量
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} {'失败':>12} {'N/A':>12} {'N/A':>12} {'N/A':>10}")
except Exception as e:
print(f"{ode_name:>15} {method_name:>15} {'错误':>12} {'N/A':>12} {'N/A':>12} {'N/A':>10}")
print()
# 可视化ODE求解结果
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]
# 绘制解析解
if ode_name == 'harmonic':
analytical_values = [harmonic_analytical(t) for t in t_eval]
ax.plot(t_eval, analytical_values, 'k--', linewidth=3, label='解析解', alpha=0.8)
else:
analytical_values = [analytical_func(t) for t in t_eval]
ax.plot(t_eval, analytical_values, 'k--', linewidth=3, label='解析解', alpha=0.8)
# 绘制数值解
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, :] # 位置分量
else:
y_values = sol.y[0, :]
ax.plot(sol.t, y_values, color=colors[i], linewidth=2,
label=f'{method_name} (误差: {method_data["error"]:.2e})', alpha=0.7)
ax.set_title(f'{description}')
ax.set_xlabel('时间 t')
ax.set_ylabel('y(t)')
ax.legend()
ax.grid(True, alpha=0.3)
plot_idx += 1
# 相空间图（仅对谐振子）
if 'harmonic' in results:
ax = axes[3]
for method_name, method_data in list(results['harmonic'].items())[:3]: # 只显示前3个方法
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)
# 解析解的相空间轨迹
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] # 简化的速度
ax.plot(y_pos_analytical, y_vel_analytical, 'k--', linewidth=3,
label='解析解', alpha=0.8)
ax.set_title('谐振子相空间轨迹')
ax.set_xlabel('位置 y')
ax.set_ylabel('速度 dy/dt')
ax.legend()
ax.grid(True, alpha=0.3)
ax.axis('equal')
plt.tight_layout()
plt.show()2. 刚性方程和高阶方程
python
# 刚性方程示例
def stiff_ode_example():
"""刚性方程示例"""
# Van der Pol振荡器 (刚性参数μ)
def van_der_pol(t, y, mu=10):
y1, y2 = y
return [y2, mu * (1 - y1**2) * y2 - y1]
# 化学反应动力学 (Robertson问题)
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振荡器'),
'robertson': (robertson, [1, 0, 0], 'Robertson化学反应')
}
# 获取刚性方程示例
stiff_examples = stiff_ode_example()
# 比较适合刚性方程的方法
stiff_methods = {
'Radau': 'Radau', # 隐式方法，适合刚性方程
'BDF': 'BDF', # 后向差分，适合刚性方程
'RK45': 'RK45' # 显式方法，对比用
}
print("刚性方程求解比较:")
print(f"{'方程':>15} {'方法':>10} {'成功':>8} {'步数':>8} {'时间(ms)':>10} {'最终时间':>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 = "是" if sol.success else "否"
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} {'否':>8} {'N/A':>8} {'N/A':>10} {'N/A':>10}")
stiff_results[eq_name][method_name] = {'success': False}
print()
# 可视化刚性方程解
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
# Van der Pol振荡器
if 'van_der_pol' in stiff_results:
# 时间序列
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振荡器 - 时间序列')
ax1.set_xlabel('时间 t')
ax1.set_ylabel('y1(t)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 相空间
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振荡器 - 相空间')
ax2.set_xlabel('y1')
ax2.set_ylabel('y2')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Robertson反应
if 'robertson' in stiff_results:
ax3 = axes[1, 0]
ax4 = axes[1, 1]
# 浓度随时间变化
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反应 - 浓度变化')
ax3.set_xlabel('时间 t')
ax3.set_ylabel('浓度 (对数尺度)')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 守恒量检查
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反应 - 守恒量偏差')
ax4.set_xlabel('时间 t')
ax4.set_ylabel('总浓度 - 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. 边值问题
python
# 边值问题求解
def boundary_value_problems():
"""边值问题示例"""
# 示例1: 二阶线性BVP
# y'' + y = 0, y(0) = 0, y(π) = 0
# 解析解: y(x) = A*sin(x), 由边界条件得 A = 0 或 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
# 示例2: 非线性BVP (Bratu问题)
# 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], '线性BVP: y\'\' + y = 0'),
'bratu': (bratu_bvp, bratu_bc, [0, 1], 'Bratu问题: y\'\' + λe^y = 0')
}
# 获取BVP示例
bvp_examples = boundary_value_problems()
print("边值问题求解:")
print(f"{'问题':>15} {'收敛':>8} {'残差':>12} {'迭代次数':>10}")
print("-" * 50)
bvp_results = {}
for bvp_name, (bvp_func, bc_func, domain, description) in bvp_examples.items():
try:
# 初始网格
x = np.linspace(domain[0], domain[1], 11)
# 初始猜测
if bvp_name == 'linear':
y_guess = np.zeros((2, x.size))
y_guess[0] = np.sin(x) # 猜测解的形状
else: # bratu
y_guess = np.zeros((2, x.size))
y_guess[0] = x * (1 - x) # 满足边界条件的初始猜测
# 求解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 = "是" if sol.success else "否"
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} {'否':>8} {'N/A':>12} {'N/A':>10}")
bvp_results[bvp_name] = {'success': False}
# 可视化BVP解
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']
# 绘制解
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='数值解')
ax.plot(sol.x, sol.y[0], 'ro', markersize=6, label='网格点')
# 添加边界条件标记
ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
ax.plot([sol.x[0], sol.x[-1]], [0, 0], 'go', markersize=8, label='边界条件')
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)
# 添加残差信息
ax.text(0.05, 0.95, f'最大残差: {result["residual"]:.2e}',
transform=ax.transAxes, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
plt.tight_layout()
plt.show()实际应用案例
1. 人口动力学模型
python
# 人口动力学建模
def population_dynamics():
"""人口动力学模型"""
# Lotka-Volterra捕食者-猎物模型
def lotka_volterra(t, y, a=1, b=0.1, c=1.5, d=0.075):
"""Lotka-Volterra方程
y[0]: 猎物数量
y[1]: 捕食者数量
"""
prey, predator = y
dprey_dt = a * prey - b * prey * predator
dpredator_dt = -c * predator + d * prey * predator
return [dprey_dt, dpredator_dt]
# 竞争模型
def competition_model(t, y, r1=1, r2=0.8, K1=100, K2=80, a12=0.5, a21=0.7):
"""种间竞争模型"""
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传染病模型
def sir_model(t, y, beta=0.3, gamma=0.1):
"""SIR传染病模型"""
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捕食模型'),
'competition': (competition_model, [50, 40], '种间竞争模型'),
'sir': (sir_model, [990, 10, 0], 'SIR传染病模型')
}
# 获取人口动力学模型
population_models = population_dynamics()
# 求解人口动力学方程
t_span = (0, 20)
t_eval = np.linspace(0, 20, 1000)
print("人口动力学模型求解:")
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}: {'成功' if sol.success else '失败'}")
except Exception as e:
print(f"{model_name:>15}: 错误 - {e}")
population_results[model_name] = {'success': False}
# 可视化人口动力学结果
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']
# 时间序列图
ax1 = axes[plot_idx]
plot_idx += 1
if model_name == 'lotka_volterra':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='猎物')
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='捕食者')
ax1.set_ylabel('数量')
elif model_name == 'competition':
ax1.plot(sol.t, sol.y[0], 'g-', linewidth=2, label='物种1')
ax1.plot(sol.t, sol.y[1], 'm-', linewidth=2, label='物种2')
ax1.set_ylabel('数量')
elif model_name == 'sir':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='易感者(S)')
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='感染者(I)')
ax1.plot(sol.t, sol.y[2], 'g-', linewidth=2, label='康复者(R)')
ax1.set_ylabel('人数')
ax1.set_title(f'{description} - 时间序列')
ax1.set_xlabel('时间')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 相空间图（除了SIR模型）
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='起点')
ax2.plot(sol.y[0, -1], sol.y[1, -1], 'ro', markersize=8, label='终点')
if model_name == 'lotka_volterra':
ax2.set_xlabel('猎物数量')
ax2.set_ylabel('捕食者数量')
ax2.set_title('Lotka-Volterra相空间轨迹')
elif model_name == 'competition':
ax2.set_xlabel('物种1数量')
ax2.set_ylabel('物种2数量')
ax2.set_title('竞争模型相空间轨迹')
ax2.legend()
ax2.grid(True, alpha=0.3)
# SIR模型的特殊分析
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']
# 计算基本再生数R0和峰值时间
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]
# 绘制感染者曲线和峰值
ax_sir.plot(sol_sir.t, sol_sir.y[1], 'r-', linewidth=3, label='感染者')
ax_sir.axvline(x=t_peak, color='orange', linestyle='--',
label=f'峰值时间: {t_peak:.1f}')
ax_sir.axhline(y=I_peak, color='orange', linestyle='--',
label=f'峰值人数: {I_peak:.0f}')
ax_sir.set_title('SIR模型 - 感染者峰值分析')
ax_sir.set_xlabel('时间')
ax_sir.set_ylabel('感染者人数')
ax_sir.legend()
ax_sir.grid(True, alpha=0.3)
# 添加统计信息
total_infected = 1000 - sol_sir.y[0, -1] # 最终易感者减少量
ax_sir.text(0.05, 0.95, f'总感染人数: {total_infected:.0f}\n感染率: {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()
# 参数敏感性分析
print("\n参数敏感性分析 (SIR模型):")
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))
# β参数影响
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('传染率β对感染者峰值的影响')
axes[0].set_xlabel('时间')
axes[0].set_ylabel('感染者人数')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# γ参数影响
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('康复率γ对感染者峰值的影响')
axes[1].set_xlabel('时间')
axes[1].set_ylabel('感染者人数')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()2. 物理系统建模
python
# 物理系统建模
def physics_systems():
"""物理系统微分方程"""
# 单摆方程
def pendulum(t, y, g=9.81, L=1.0, damping=0.1):
"""阻尼单摆方程
θ'' + (damping/m)θ' + (g/L)sin(θ) = 0
"""
theta, omega = y
return [omega, -damping * omega - (g/L) * np.sin(theta)]
# 弹簧-质量-阻尼系统
def spring_mass_damper(t, y, m=1.0, k=4.0, c=0.5):
"""弹簧-质量-阻尼系统
mx'' + cx' + kx = 0
"""
x, v = y
return [v, -(c/m) * v - (k/m) * x]
# RLC电路
def rlc_circuit(t, y, R=1.0, L=1.0, C=1.0, V0=0):
"""RLC电路方程
L(dI/dt) + R*I + Q/C = V0
dQ/dt = I
"""
Q, I = y # 电荷和电流
return [I, (V0 - R*I - Q/C) / L]
# 洛伦兹方程（混沌系统）
def lorenz(t, y, sigma=10, rho=28, beta=8/3):
"""洛伦兹方程"""
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], '阻尼单摆'),
'spring_mass': (spring_mass_damper, [1, 0], '弹簧-质量-阻尼系统'),
'rlc_circuit': (rlc_circuit, [0, 0], 'RLC电路'),
'lorenz': (lorenz, [1, 1, 1], '洛伦兹混沌系统')
}
# 获取物理系统模型
physics_models = physics_systems()
# 求解物理系统方程
t_span_physics = (0, 10)
t_eval_physics = np.linspace(0, 10, 2000)
print("物理系统建模:")
physics_results = {}
for model_name, (model_func, y0, description) in physics_models.items():
try:
if model_name == 'lorenz':
# 洛伦兹系统需要更长的时间
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}: {'成功' if sol.success else '失败'}")
except Exception as e:
print(f"{model_name:>15}: 错误 - {e}")
physics_results[model_name] = {'success': False}
# 可视化物理系统结果
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']
# 时间序列图
ax1 = axes[plot_idx]
plot_idx += 1
if model_name == 'pendulum':
ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='角度θ')
ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='角速度ω')
ax1.set_ylabel('角度/角速度')
elif model_name == 'spring_mass':
ax1.plot(sol.t, sol.y[0], 'g-', linewidth=2, label='位移x')
ax1.plot(sol.t, sol.y[1], 'm-', linewidth=2, label='速度v')
ax1.set_ylabel('位移/速度')
elif model_name == 'rlc_circuit':
ax1.plot(sol.t, sol.y[0], 'c-', linewidth=2, label='电荷Q')
ax1.plot(sol.t, sol.y[1], 'orange', linewidth=2, label='电流I')
ax1.set_ylabel('电荷/电流')
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('坐标值')
ax1.set_title(f'{description} - 时间序列')
ax1.set_xlabel('时间')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 相空间图
if plot_idx < 8:
ax2 = axes[plot_idx]
plot_idx += 1
if model_name == 'lorenz':
# 3D轨迹投影到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('洛伦兹吸引子 (x-y投影)')
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='起点')
if model_name == 'pendulum':
ax2.set_xlabel('角度θ')
ax2.set_ylabel('角速度ω')
ax2.set_title('单摆相空间轨迹')
elif model_name == 'spring_mass':
ax2.set_xlabel('位移x')
ax2.set_ylabel('速度v')
ax2.set_title('弹簧系统相空间轨迹')
elif model_name == 'rlc_circuit':
ax2.set_xlabel('电荷Q')
ax2.set_ylabel('电流I')
ax2.set_title('RLC电路相空间轨迹')
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 洛伦兹系统的3D可视化
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轨迹
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('洛伦兹吸引子 3D轨迹')
# x-z投影
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('洛伦兹吸引子 (x-z投影)')
ax2.grid(True, alpha=0.3)
# y-z投影
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('洛伦兹吸引子 (y-z投影)')
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
## 性能优化技巧
### 1. 算法选择和参数调优
```python
# 性能优化示例
def performance_optimization():
"""积分和ODE求解的性能优化"""
# 测试函数
def test_integration_performance():
"""积分性能测试"""
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("积分方法性能比较:")
print(f"{'函数':>12} {'方法':>15} {'结果':>12} {'时间(ms)':>10} {'精度':>12}")
print("-" * 70)
for func_name, func in functions.items():
# 参考结果
ref_result, _ = quad(func, 0, 1)
for method_name, method in methods.items():
times = []
for _ in range(10): # 多次测试取平均
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求解性能测试"""
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求解方法性能比较:")
print(f"{'方程类型':>12} {'方法':>10} {'成功':>6} {'步数':>8} {'时间(ms)':>10} {'精度':>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 = "是" if sol.success else "否"
n_steps = len(sol.t)
# 简单的精度估计（与最精确方法比较）
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} {'否':>6} {'N/A':>8} {'N/A':>10} {'N/A':>12}")
print()
# 执行性能测试
test_integration_performance()
test_ode_performance()
# 执行性能优化测试
performance_optimization()2. 内存优化和并行计算
python
# 内存优化和并行计算
def memory_and_parallel_optimization():
"""内存优化和并行计算示例"""
# 大规模积分的内存优化
def memory_efficient_integration():
"""内存高效的积分计算"""
def large_array_function(x):
# 模拟需要大量内存的函数
return np.sum(np.sin(np.arange(1000) * x))
# 分块积分
def chunked_integration(func, a, b, n_chunks=10):
"""分块积分以节省内存"""
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
# 比较内存使用
print("内存优化积分比较:")
# 直接积分
start_time = time.time()
direct_result, _ = quad(large_array_function, 0, 1)
direct_time = time.time() - start_time
# 分块积分
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_result:.6f}, 时间: {direct_time*1000:.2f}ms")
print(f"分块积分: {chunked_result:.6f}, 时间: {chunked_time*1000:.2f}ms")
print(f"相对误差: {abs(direct_result - chunked_result)/abs(direct_result):.2e}")
print()
# 并行ODE求解
def parallel_ode_solving():
"""并行ODE求解示例"""
def parameter_study_ode(t, y, param):
"""参数化的ODE"""
return [-param * y[0], y[0] - y[1]]
# 参数范围
parameters = np.linspace(0.1, 10, 20)
# 串行求解
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_time*1000:.2f}ms")
print(f"成功求解: {np.sum(~np.isnan(serial_results))}/{len(parameters)}")
return serial_results, parameters
# 执行优化测试
memory_efficient_integration()
results, params = parallel_ode_solving()
# 可视化参数研究结果
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(params, results, 'bo-', linewidth=2, markersize=6)
plt.xlabel('参数值')
plt.ylabel('最终状态 y[0](T)')
plt.title('参数敏感性分析')
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
# 显示收敛性
valid_results = np.array(results)[~np.isnan(results)]
valid_params = params[~np.isnan(results)]
if len(valid_results) > 1:
# 计算梯度
gradient = np.gradient(valid_results, valid_params)
plt.plot(valid_params, gradient, 'ro-', linewidth=2, markersize=6)
plt.xlabel('参数值')
plt.ylabel('敏感性 dy/dp')
plt.title('参数敏感性梯度')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 执行内存和并行优化
memory_and_parallel_optimization()
## 小结
SciPy的积分和微分方程模块提供了强大而灵活的数值计算工具：
### 积分功能
- **一维积分**: `quad`、`simpson`、`trapz`等方法适用于不同精度需求
- **多重积分**: `dblquad`、`tplquad`、`nquad`支持复杂区域积分
- **特殊积分**: 处理奇异点、无穷区间、振荡函数等特殊情况
### 微分方程求解
- **初值问题**: `solve_ivp`提供多种高精度求解器
- **边值问题**: `solve_bvp`处理复杂边界条件
- **刚性方程**: `Radau`、`BDF`等隐式方法处理刚性系统
### 应用领域
- **科学计算**: 物理建模、工程分析、数值仿真
- **生物数学**: 人口动力学、传染病模型、生态系统
- **工程应用**: 控制系统、信号处理、电路分析
### 性能优化
- **算法选择**: 根据问题特性选择合适的求解器
- **参数调优**: 合理设置容差、步长等参数
- **内存管理**: 分块处理大规模问题
- **并行计算**: 参数研究和蒙特卡罗模拟
掌握这些工具和技巧，可以高效解决各种数值积分和微分方程问题，为科学研究和工程应用提供强有力的计算支持。
## 练习题
1. **积分计算**: 计算 ∫₀^π sin(x)cos(x)dx，比较解析解和数值解的精度
2. **二重积分**: 计算圆形区域 x²+y²≤1 上函数 f(x,y)=e^(-(x²+y²)) 的积分
3. **ODE求解**: 求解 Lotka-Volterra 方程，分析不同初值对解的影响
4. **刚性方程**: 求解 Robertson 化学反应方程，比较不同求解器的性能
5. **边值问题**: 求解二阶边值问题 y''+λy=0, y(0)=y(π)=0，找出特征值λ
6. **参数研究**: 对 SIR 模型进行参数敏感性分析，研究 β 和 γ 对疫情峰值的影响
7. **混沌系统**: 分析洛伦兹方程的混沌行为，绘制不同参数下的吸引子
8. **性能优化**: 实现一个自适应积分算法，并与 scipy.integrate.quad 比较性能