python 求解二维随机半线性波动方程

在Python中，你可以使用`Scipy`的`solve_ivp`函数结合`Euler-Maruyama`或者`Midpoint`等方法来求解二维的随机半线性波动方程。这个方程通常形式如下：

\[ \frac{dU}{dt} = a(U, V) + b(U, V)dW_t \]
\[ \frac{dV}{dt} = c(U, V) \]

其中 \( U \) 和 \( V \) 是状态变量，\( a(U, V), b(U, V), c(U, V) \) 是函数，\( dW_t \) 表示布朗运动噪声。

下面是一个基于`Euler-Maruyama`方法的基本示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 假设我们有简单的二维半线性方程，例如：
def sde(t, y, dt, random_var):  
    U, V = y
    dU = -V * dt + np.sqrt(U) * random_var  # 二次增长项
    dV = U * dt  # 线性项
    return [dU, dV]

# 设置参数
a = 1.0
b = 0.5  # 白噪声系数
sigma = 0.1  # 噪声的标准差
t0, tf = 0., 1.  # 时间范围
y0 = [1., 0.]  # 初始值
dt = 0.01  # 步长

# 创建随机噪声序列
random_sequence = np.random.normal(loc=0, scale=sigma, size=int((tf - t0) // dt))

# 使用Euler-Maruyama方法求解
sol = solve_ivp(sde, (t0, tf), y0, args=(dt, random_sequence), method='Euler-Maruyama', dense_output=True)

# 绘制解
t = sol.t
U = sol.y[0]
V = sol.y[1]
plt.plot(t, U, label='U')
plt.plot(t, V, label='V')
plt.xlabel('time')
plt.ylabel('solution')
plt.legend()