当前位置： 首页 > news >正文

Python实例题：Python计算常微分方程

news 来源：原创 2025/6/14 7:13:49

目录

Python实例题

题目

处理复杂概率分布的方法

自定义概率分布：使用scipy.stats创建自定义分布

使用混合分布：将多个简单分布组合成复杂分布

接受 - 拒绝采样：从复杂分布中生成样本

MCMC 方法：使用马尔可夫链蒙特卡洛采样

Python 计算常微分方程

使用 scipy.integrate.odeint：最常用的 ODE 求解器

使用 scipy.integrate.solve_ivp：更现代的 ODE 求解器

求解系统 ODE：多个相互关联的微分方程

高阶 ODE 转换为一阶系统：处理二阶及以上的 ODE

边界值问题 (BVP)：使用 scipy.integrate.solve_bvp

处理复杂概率分布与 ODE 结合的随机微分方程

Python实例题

题目

Python计算常微分方程

处理复杂概率分布的方法

在计算随机过程时，处理复杂概率分布是常见的挑战。以下是几种有效的方法：

  • 自定义概率分布：使用scipy.stats创建自定义分布

from scipy import stats
import numpy as np
import matplotlib.pyplot as plt# 定义自定义概率密度函数
def custom_pdf(x):return np.exp(-x**2/2) * (1 + 0.5 * np.sin(2*x)) / np.sqrt(2*np.pi)# 创建自定义分布
class CustomDistribution(stats.rv_continuous):def _pdf(self, x):return custom_pdf(x)# 实例化分布
custom_dist = CustomDistribution(name='custom')# 生成随机样本
samples = custom_dist.rvs(size=1000)# 可视化
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, label='样本分布')
x = np.linspace(-4, 4, 1000)
plt.plot(x, custom_pdf(x), 'r-', label='理论PDF')
plt.legend()
plt.title('自定义概率分布')
plt.show()

  • 使用混合分布：将多个简单分布组合成复杂分布

# 混合高斯分布
def mixture_gaussian(x, weights, means, stds):result = 0for w, mu, sigma in zip(weights, means, stds):result += w * stats.norm.pdf(x, mu, sigma)return result# 参数
weights = [0.3, 0.5, 0.2]
means = [-2, 0, 3]
stds = [0.8, 1.2, 0.7]# 生成随机样本
n_samples = 1000
component = np.random.choice(len(weights), size=n_samples, p=weights)
samples = np.zeros(n_samples)for i in range(n_samples):samples[i] = np.random.normal(means[component[i]], stds[component[i]])# 可视化
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, label='样本分布')
x = np.linspace(-6, 6, 1000)
plt.plot(x, mixture_gaussian(x, weights, means, stds), 'r-', label='理论PDF')
plt.legend()
plt.title('混合高斯分布')
plt.show()

  • 接受 - 拒绝采样：从复杂分布中生成样本

def complex_pdf(x):# 复杂分布的PDFreturn 0.3 * np.exp(-0.2 * x**2) + 0.7 * np.exp(-0.2 * (x-10)**2)def sampling_complex(n_samples, proposal_dist=stats.norm, c=1.5):samples = []while len(samples) < n_samples:# 从提议分布中采样x = proposal_dist.rvs()# 计算接受概率accept_prob = complex_pdf(x) / (c * proposal_dist.pdf(x))# 决定是否接受if np.random.random() < accept_prob:samples.append(x)return np.array(samples)# 生成样本
samples = sampling_complex(1000)# 可视化
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, label='样本分布')
x = np.linspace(-5, 15, 1000)
plt.plot(x, complex_pdf(x), 'r-', label='理论PDF')
plt.legend()
plt.title('接受-拒绝采样')
plt.show()

  • MCMC 方法：使用马尔可夫链蒙特卡洛采样

def metropolis_hastings(target_pdf, n_samples, proposal_std=1.0, x0=0):samples = [x0]current_x = x0for i in range(n_samples - 1):# 提议新样本proposed_x = current_x + np.random.normal(0, proposal_std)# 计算接受概率ratio = target_pdf(proposed_x) / target_pdf(current_x)accept_prob = min(1, ratio)# 决定是否接受if np.random.random() < accept_prob:current_x = proposed_xsamples.append(current_x)return np.array(samples)# 使用Metropolis-Hastings从复杂分布采样
samples = metropolis_hastings(complex_pdf, 5000)
samples = samples[1000:]  # 丢弃burn-in期# 可视化
plt.figure(figsize=(10, 6))
plt.hist(samples, bins=50, density=True, alpha=0.6, label='样本分布')
x = np.linspace(-5, 15, 1000)
plt.plot(x, complex_pdf(x), 'r-', label='理论PDF')
plt.legend()
plt.title('MCMC采样')
plt.show()

Python 计算常微分方程

Python 提供了多种方法来求解常微分方程 (ODEs)。下面介绍主要的方法和库：

  • 使用 scipy.integrate.odeint：最常用的 ODE 求解器

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt# 定义微分方程 dy/dt = f(y, t)
def model(y, t, k):dydt = -k * yreturn dydt# 初始条件
y0 = 5# 时间点
t = np.linspace(0, 20, 100)# 求解ODE
k = 0.1
y = odeint(model, y0, t, args=(k,))# 可视化结果
plt.figure(figsize=(10, 6))
plt.plot(t, y, 'b-', linewidth=2)
plt.xlabel('时间')
plt.ylabel('y(t)')
plt.title('一阶常微分方程 dy/dt = -ky 的解')
plt.grid(True)
plt.show()

  • 使用 scipy.integrate.solve_ivp：更现代的 ODE 求解器

from scipy.integrate import solve_ivp# 定义微分方程 dy/dt = f(t, y)
def model(t, y, k):return -k * y# 时间跨度
t_span = (0, 20)
t_eval = np.linspace(0, 20, 100)# 初始条件
y0 = [5]# 求解ODE
k = 0.1
sol = solve_ivp(lambda t, y: model(t, y, k), t_span, y0, t_eval=t_eval)# 可视化结果
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], 'b-', linewidth=2)
plt.xlabel('时间')
plt.ylabel('y(t)')
plt.title('使用solve_ivp求解一阶常微分方程')
plt.grid(True)
plt.show()

  • 求解系统 ODE：多个相互关联的微分方程

# 定义Lotka-Volterra捕食者-猎物模型
def lotka_volterra(t, z, a, b, c, d):x, y = zdxdt = a * x - b * x * ydydt = -c * y + d * x * yreturn [dxdt, dydt]# 参数
a = 1.0
b = 0.1
c = 1.0
d = 0.02# 初始条件
z0 = [10, 5]  # 初始猎物数量和捕食者数量# 时间跨度
t_span = (0, 15)
t_eval = np.linspace(0, 15, 1000)# 求解ODE
sol = solve_ivp(lambda t, z: lotka_volterra(t, z, a, b, c, d), t_span, z0, t_eval=t_eval, method='RK45')# 可视化结果
plt.figure(figsize=(12, 6))
plt.plot(sol.t, sol.y[0], 'b-', label='猎物数量')
plt.plot(sol.t, sol.y[1], 'r-', label='捕食者数量')
plt.xlabel('时间')
plt.ylabel('种群数量')
plt.title('Lotka-Volterra捕食者-猎物模型')
plt.legend()
plt.grid(True)
plt.show()

  • 高阶 ODE 转换为一阶系统：处理二阶及以上的 ODE

# 阻尼振荡器方程: y'' + 2*zeta*omega_n*y' + omega_n^2*y = 0
def oscillator(t, z, omega_n, zeta):y, v = z  # y是位置，v是速度dydt = vdvdt = -omega_n**2 * y - 2 * zeta * omega_n * vreturn [dydt, dvdt]# 参数
omega_n = 2.0  # 自然频率
zeta = 0.1     # 阻尼比# 初始条件
z0 = [1.0, 0.0]  # 初始位置和速度# 时间跨度
t_span = (0, 10)
t_eval = np.linspace(0, 10, 1000)# 求解ODE
sol = solve_ivp(lambda t, z: oscillator(t, z, omega_n, zeta), t_span, z0, t_eval=t_eval)# 可视化结果
plt.figure(figsize=(12, 6))
plt.plot(sol.t, sol.y[0], 'b-', label='位置')
plt.plot(sol.t, sol.y[1], 'r-', label='速度')
plt.xlabel('时间')
plt.ylabel('状态')
plt.title('阻尼振荡器')
plt.legend()
plt.grid(True)
plt.show()

  • 边界值问题 (BVP)：使用 scipy.integrate.solve_bvp

from scipy.integrate import solve_bvp# 定义边界值问题 y'' + y = 0
def bvp_ode(x, y):return np.vstack((y[1], -y[0]))# 定义边界条件 y(0) = 0, y(pi/2) = 1
def bvp_bc(ya, yb):return np.array([ya[0], yb[0] - 1])# 初始网格和猜测
x = np.linspace(0, np.pi/2, 5)
y = np.zeros((2, x.size))
y[0] = x  # 初始猜测# 求解BVP
sol = solve_bvp(bvp_ode, bvp_bc, x, y)# 可视化结果
x_plot = np.linspace(0, np.pi/2, 100)
y_plot = sol.sol(x_plot)[0]plt.figure(figsize=(10, 6))
plt.plot(x_plot, y_plot, 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title('边界值问题 y\'\' + y = 0 的解')
plt.grid(True)
plt.show()

处理复杂概率分布与 ODE 结合的随机微分方程

当需要求解随机微分方程 (SDE) 时，上述技术可以结合使用：

# 几何布朗运动的随机微分方程 dS = mu*S*dt + sigma*S*dW
def geometric_brownian_motion(S0, mu, sigma, T, N, M):"""模拟几何布朗运动参数:S0: 初始价格mu: 漂移率sigma: 波动率T: 总时间N: 时间步数M: 模拟路径数"""dt = T / Nt = np.linspace(0, T, N+1)# 生成随机增量dW = np.random.normal(0, np.sqrt(dt), (M, N))# 计算价格路径S = np.zeros((M, N+1))S[:, 0] = S0for i in range(N):# Euler-Maruyama方法S[:, i+1] = S[:, i] * (1 + mu * dt + sigma * dW[:, i])return t, S# 参数
S0 = 100      # 初始价格
mu = 0.05     # 漂移率
sigma = 0.2   # 波动率
T = 1.0       # 总时间
N = 252       # 交易日数
M = 10        # 模拟路径数# 模拟
t, S = geometric_brownian_motion(S0, mu, sigma, T, N, M)# 可视化
plt.figure(figsize=(12, 6))
for i in range(M):plt.plot(t, S[i], label=f'路径 {i+1}')plt.xlabel('时间')
plt.ylabel('价格')
plt.title('几何布朗运动模拟')
plt.legend()
plt.grid(True)
plt.show()

相关文章：

  • 电路图识图基础知识-电动机的保护电路保护方式（二十六）
  • 胶囊网络破解图像旋转不变性难题 ——从空间关系到姿态矩阵的几何深度学习革命
  • C++进阶—C++中的继承
  • BIO网络通信基础（TCP协议）
  • 基于 STL+VMD 二次分解的 Informer-LSTM 并行预测模型详解与案例
  • USB串口通信、握手协议、深度学习等技术要点
  • L1-078 吉老师的回归
  • vue3项目使用静态资源
  • 以楼宇自控关键技术，夯实现代低碳建筑发展重要基础
  • linux命令 root
  • PGCP：用于比较基因组学的植物基因组综合数据库-文献精读144
  • colmap编译到本地【Ubuntu22.04】
  • HarmonyOS5 运动健康app(二)：健康跑步（附代码）
  • Text2SQL、Text2API基础
  • Java实现Pdf转Word
  • Excel大厂自动化报表实战(互联网金融-数据分析周报制作中)
  • 逆向入门（3）程序逆向篇-Acid Bytes.2
  • ＜teleport＞ 是 Vue 3 引入的一个内置组件，用于在 DOM 中移动组件的渲染位置，但保持组件的逻辑作用域不变
  • 【62 Pandas+Pyecharts | 智联招聘大数据岗位数据分析可视化】
  • 【舞蹈】FineDance
  • 河北省建设厅官方网站 官网/公司运营策划方案
  • 做商城网站需要备案什么域名/手机怎么建网站
  • 网站建设明薇通网络/做一个微信小程序需要多少钱
  • 如何用github做网站/黑龙江网络推广好做吗
  • 网站建设单页面推广模板/营销型网站建设套餐
  • 中国建设招标网是什么网站/网站关键词提升