# Python常微分方程 - 完整代码示例# 常微分方程描述系统随时间的演化广泛应用于物理、生物和经济建模import numpy as npimport matplotlib.pyplot as pltfrom scipy import integrate# 1. 单变量 ODE指数衰减模型# 方程: dy/dt -ky, y(0) y0, 解析解: y(t) y0 * exp(-kt)def exponential_decay(t, y, k0.5):指数衰减的右侧函数 dy/dt -k*yreturn -k * y# 求解 ODEt_span (0, 10) # 时间区间y0 [5.0] # 初始条件# 使用 RK45 方法求解sol integrate.solve_ivp(exponential_decay, t_span, y0,methodRK45, dense_outputTrue,max_step0.1 # 最大步长控制精度)# 在均匀时间点上评估解t_eval np.linspace(0, 10, 100)y_sol sol.sol(t_eval)[0] # 使用 dense_output 插值y_analytic 5.0 * np.exp(-0.5 * t_eval) # 解析解print(单变量 ODE 求解:)print(ft5 处的数值解: {sol.sol(5)[0]:.6f})print(ft5 处的解析解: {5*np.exp(-0.5*5):.6f})print(f最大误差: {np.max(np.abs(y_sol - y_analytic)):.2e})# 2. 一阶 ODE 系统Lotka-Volterra 捕食者-猎物模型# 方程: dx/dt ax - bxy, dy/dt -cy dxy# x: 猎物数量, y: 捕食者数量def lotka_volterra(t, z, a, b, c, d):Lotka-Volterra 捕食者-猎物系统x, y z # z [x, y] 是状态向量dxdt a * x - b * x * y # 猎物增长率dydt -c * y d * x * y # 捕食者增长率return [dxdt, dydt]# 参数和初始条件a, b, c, d 1.5, 1.0, 3.0, 1.0 # 模型参数z0 [10.0, 5.0] # 初始种群数量 [猎物, 捕食者]# 求解系统sol_lv integrate.solve_ivp(lotka_volterra, (0, 20), z0,args(a, b, c, d),methodRK45, dense_outputTrue,max_step0.05)t_lv np.linspace(0, 20, 500)x_lv, y_lv sol_lv.sol(t_lv)print(f\nLotka-Volterra 系统:)print(f猎物数量范围: [{np.min(x_lv):.1f}, {np.max(x_lv):.1f}])print(f捕食者数量范围: [{np.min(y_lv):.1f}, {np.max(y_lv):.1f}])# 3. Lorenz 系统混沌系统# 方程: dx/dt σ(y-x), dy/dt x(ρ-z)-y, dz/dt xy - βzdef lorenz_system(t, state, sigma, rho, beta):Lorenz 混沌系统x, y, z statedx sigma * (y - x)dy x * (rho - z) - ydz x * y - beta * zreturn [dx, dy, dz]# 经典参数混沌吸引子sigma, rho, beta 10, 28, 8/3state0 [1.0, 1.0, 1.0]# 使用 BDF 方法刚性系统适用sol_lorenz integrate.solve_ivp(lorenz_system, (0, 40), state0,args(sigma, rho, beta),methodBDF, dense_outputTrue,rtol1e-9, atol1e-12 # 高精度要求)t_lor np.linspace(0, 40, 5000)x_lor, y_lor, z_lor sol_lorenz.sol(t_lor)print(f\nLorenz 系统混沌:)print(fx 范围: [{np.min(x_lor):.1f}, {np.max(x_lor):.1f}])print(fy 范围: [{np.min(y_lor):.1f}, {np.max(y_lor):.1f}])print(fz 范围: [{np.min(z_lor):.1f}, {np.max(z_lor):.1f}])# 4. 刚性 ODE范德波尔振荡器# 当 μ 1 时系统变刚需使用 BDF 或 Radau 方法def van_der_pol(t, state, mu):范德波尔振荡器 y - μ(1-y²)y y 0y, v state # v dy/dtdy vdv mu * (1 - y**2) * v - yreturn [dy, dv]mu 1000 # 大 μ 使系统变刚sol_vdp integrate.solve_ivp(van_der_pol, (0, 3000), [2, 0],args(mu,), methodRadau,rtol1e-6, atol1e-8)print(f\n范德波尔振荡器 (μ{mu}):)print(f求解点数: {len(sol_vdp.t)})# 5. 方法对比不同求解器的性能import timedef compare_methods(ode_func, t_span, y0, methods):对比不同 ODE 求解器的性能for method in methods:start time.time()sol integrate.solve_ivp(ode_func, t_span, y0, methodmethod,rtol1e-6, atol1e-8)elapsed time.time() - startprint(f {method:12s}: {len(sol.t):6d} 步, {elapsed:.4f}s, 状态{sol.status})print(f\n求解器性能对比 (Lotka-Volterra):)compare_methods(lambda t, z: lotka_volterra(t, z, a, b, c, d),(0, 20), z0,[RK45, RK23, DOP853, Radau, BDF])# 6. 事件检测使用 events 参数# 检测种群数量何时达到特定值def event_cross_x(t, z):事件函数当 x 15 时触发x, y zreturn x - 15event_cross_x.terminal True # 触发后停止积分event_cross_x.direction 1 # 只检测正向穿越sol_event integrate.solve_ivp(lambda t, z: lotka_volterra(t, z, a, b, c, d),(0, 20), z0, eventsevent_cross_x,max_step0.1)if sol_event.t_events[0].size 0:print(f\n事件检测: x15 发生在 t{sol_event.t_events[0][0]:.4f})# 7. 相图分析和轨道稳定性# 计算多个不同初始条件的轨道print(f\n不同初始条件的轨道终点:)for x0_init in [5, 10, 15]:sol_test integrate.solve_ivp(lambda t, z: lotka_volterra(t, z, a, b, c, d),(0, 50), [x0_init, 5.0], methodRK45,max_step0.1)print(f x0{x0_init:3d}: 终点 x{sol_test.y[0,-1]:.2f}, y{sol_test.y[1,-1]:.2f})print(\n常微分方程总结选对求解器和方法至关重要)print(RK45 通用BDF/Radau 适合刚性系统事件检测增强实用性)
Python常微分方程
# Python常微分方程 - 完整代码示例# 常微分方程描述系统随时间的演化广泛应用于物理、生物和经济建模import numpy as npimport matplotlib.pyplot as pltfrom scipy import integrate# 1. 单变量 ODE指数衰减模型# 方程: dy/dt -ky, y(0) y0, 解析解: y(t) y0 * exp(-kt)def exponential_decay(t, y, k0.5):指数衰减的右侧函数 dy/dt -k*yreturn -k * y# 求解 ODEt_span (0, 10) # 时间区间y0 [5.0] # 初始条件# 使用 RK45 方法求解sol integrate.solve_ivp(exponential_decay, t_span, y0,methodRK45, dense_outputTrue,max_step0.1 # 最大步长控制精度)# 在均匀时间点上评估解t_eval np.linspace(0, 10, 100)y_sol sol.sol(t_eval)[0] # 使用 dense_output 插值y_analytic 5.0 * np.exp(-0.5 * t_eval) # 解析解print(单变量 ODE 求解:)print(ft5 处的数值解: {sol.sol(5)[0]:.6f})print(ft5 处的解析解: {5*np.exp(-0.5*5):.6f})print(f最大误差: {np.max(np.abs(y_sol - y_analytic)):.2e})# 2. 一阶 ODE 系统Lotka-Volterra 捕食者-猎物模型# 方程: dx/dt ax - bxy, dy/dt -cy dxy# x: 猎物数量, y: 捕食者数量def lotka_volterra(t, z, a, b, c, d):Lotka-Volterra 捕食者-猎物系统x, y z # z [x, y] 是状态向量dxdt a * x - b * x * y # 猎物增长率dydt -c * y d * x * y # 捕食者增长率return [dxdt, dydt]# 参数和初始条件a, b, c, d 1.5, 1.0, 3.0, 1.0 # 模型参数z0 [10.0, 5.0] # 初始种群数量 [猎物, 捕食者]# 求解系统sol_lv integrate.solve_ivp(lotka_volterra, (0, 20), z0,args(a, b, c, d),methodRK45, dense_outputTrue,max_step0.05)t_lv np.linspace(0, 20, 500)x_lv, y_lv sol_lv.sol(t_lv)print(f\nLotka-Volterra 系统:)print(f猎物数量范围: [{np.min(x_lv):.1f}, {np.max(x_lv):.1f}])print(f捕食者数量范围: [{np.min(y_lv):.1f}, {np.max(y_lv):.1f}])# 3. Lorenz 系统混沌系统# 方程: dx/dt σ(y-x), dy/dt x(ρ-z)-y, dz/dt xy - βzdef lorenz_system(t, state, sigma, rho, beta):Lorenz 混沌系统x, y, z statedx sigma * (y - x)dy x * (rho - z) - ydz x * y - beta * zreturn [dx, dy, dz]# 经典参数混沌吸引子sigma, rho, beta 10, 28, 8/3state0 [1.0, 1.0, 1.0]# 使用 BDF 方法刚性系统适用sol_lorenz integrate.solve_ivp(lorenz_system, (0, 40), state0,args(sigma, rho, beta),methodBDF, dense_outputTrue,rtol1e-9, atol1e-12 # 高精度要求)t_lor np.linspace(0, 40, 5000)x_lor, y_lor, z_lor sol_lorenz.sol(t_lor)print(f\nLorenz 系统混沌:)print(fx 范围: [{np.min(x_lor):.1f}, {np.max(x_lor):.1f}])print(fy 范围: [{np.min(y_lor):.1f}, {np.max(y_lor):.1f}])print(fz 范围: [{np.min(z_lor):.1f}, {np.max(z_lor):.1f}])# 4. 刚性 ODE范德波尔振荡器# 当 μ 1 时系统变刚需使用 BDF 或 Radau 方法def van_der_pol(t, state, mu):范德波尔振荡器 y - μ(1-y²)y y 0y, v state # v dy/dtdy vdv mu * (1 - y**2) * v - yreturn [dy, dv]mu 1000 # 大 μ 使系统变刚sol_vdp integrate.solve_ivp(van_der_pol, (0, 3000), [2, 0],args(mu,), methodRadau,rtol1e-6, atol1e-8)print(f\n范德波尔振荡器 (μ{mu}):)print(f求解点数: {len(sol_vdp.t)})# 5. 方法对比不同求解器的性能import timedef compare_methods(ode_func, t_span, y0, methods):对比不同 ODE 求解器的性能for method in methods:start time.time()sol integrate.solve_ivp(ode_func, t_span, y0, methodmethod,rtol1e-6, atol1e-8)elapsed time.time() - startprint(f {method:12s}: {len(sol.t):6d} 步, {elapsed:.4f}s, 状态{sol.status})print(f\n求解器性能对比 (Lotka-Volterra):)compare_methods(lambda t, z: lotka_volterra(t, z, a, b, c, d),(0, 20), z0,[RK45, RK23, DOP853, Radau, BDF])# 6. 事件检测使用 events 参数# 检测种群数量何时达到特定值def event_cross_x(t, z):事件函数当 x 15 时触发x, y zreturn x - 15event_cross_x.terminal True # 触发后停止积分event_cross_x.direction 1 # 只检测正向穿越sol_event integrate.solve_ivp(lambda t, z: lotka_volterra(t, z, a, b, c, d),(0, 20), z0, eventsevent_cross_x,max_step0.1)if sol_event.t_events[0].size 0:print(f\n事件检测: x15 发生在 t{sol_event.t_events[0][0]:.4f})# 7. 相图分析和轨道稳定性# 计算多个不同初始条件的轨道print(f\n不同初始条件的轨道终点:)for x0_init in [5, 10, 15]:sol_test integrate.solve_ivp(lambda t, z: lotka_volterra(t, z, a, b, c, d),(0, 50), [x0_init, 5.0], methodRK45,max_step0.1)print(f x0{x0_init:3d}: 终点 x{sol_test.y[0,-1]:.2f}, y{sol_test.y[1,-1]:.2f})print(\n常微分方程总结选对求解器和方法至关重要)print(RK45 通用BDF/Radau 适合刚性系统事件检测增强实用性)
