带跳的随机微分方程的代码
时间: 2023-05-20 20:06:58 浏览: 98
以下是一个带跳的随机微分方程的代码示例:
```python
import numpy as np
import matplotlib.pyplot as plt
# 定义微分方程的参数
mu = 0.1
sigma = 0.2
lam = 0.5
delta = 0.1
# 定义初始值和时间步长
x0 = 1
t0 = 0
dt = 0.01
T = 10
# 定义随机跳跃的时间步长
dt_jump = 0.1
# 定义随机跳跃的概率分布
p = np.array([0.5, 0.5])
s = np.array([-1, 1])
# 定义微分方程的函数
def f(x):
return mu * x
def g(x):
return sigma * x
# 定义随机跳跃的函数
def jump(x):
return lam * x
# 初始化数组
n = int(T / dt)
n_jump = int(T / dt_jump)
t = np.linspace(t0, T, n+1)
x = np.zeros(n+1)
x_jump = np.zeros(n_jump+1)
# 生成随机数
r = np.random.normal(size=n)
r_jump = np.random.choice(s, size=n_jump, p=p)
# 循环计算微分方程的值
x[0] = x0
for i in range(n):
x[i+1] = x[i] + f(x[i]) * dt + g(x[i]) * np.sqrt(dt) * r[i]
if (i+1) % (dt_jump / dt) == 0:
x_jump[(i+1) // (dt_jump / dt)] = x_jump[(i+1) // (dt_jump / dt) - 1] + jump(x_jump[(i+1) // (dt_jump / dt) - 1]) * dt_jump + delta * r_jump[(i+1) // (dt_jump / dt) - 1]
else:
x_jump[(i+1) // (dt_jump / dt)] = x_jump[(i+1) // (dt_jump / dt) - 1]
# 绘制图像
plt.plot(t, x, label='diffusion')
plt.plot(np.linspace(t0, T, n_jump+1), x_jump, label='jump')
plt.legend()
plt.show()
```
这个代码使用了欧拉方法来数值求解带跳的随机微分方程,其中随机跳跃的部分使用了泊松过程来模拟。