🎯要点
- 基于随机动态行为受化学主方程控制，定量分析单细胞瞬态效应。
- 确定性常微分方程描述双稳态和滞后现象。
- 通过随机性偏微分方程描述出暂时性滞后会逐渐达到平稳状态，并利用熵方法或截断方法计算平衡收敛速度的估计值。
- 随机定量分析模型使用最小二乘法从时间依赖性分布中找到最佳参数集，在执行模拟以估计收敛速度。
🍪语言内容分比
🍇MATLAB确定性常微分方程和随机偏微分方程
在 MATLAB 中，确定性常微分方程（ODE）可以通过多种数值方法进行求解，最常用的是 MATLAB 提供的 ode45、ode23 等函数。这些函数使用不同的数值算法来求解初值问题，并且适用于不同类型的常微分方程。
下面，我们介绍一些常用方法来解确定性常微分方程。
1. 基本求解器：ode45
ode45 是 MATLAB 中最常用的 ODE 求解器，适用于大多数非刚性 ODE。ode45 基于 Runge-Kutta 方法，对函数进行迭代求解。
下面是一个例子，求解简单的常微分方程 ( \frac{dy}{dt} = -2y )，并且初始条件为 ( y(0) = 1 )。
% 定义微分方程
dydt = @(t, y) -2 * y;
% 设置时间范围
tspan = [0 5]; % 从0到5秒
% 设置初始条件
y0 = 1;
% 使用 ode45 求解 ODE
[t, y] = ode45(dydt, tspan, y0);
% 绘图
plot(t, y);
title('解 \frac{dy}{dt} = -2y');
xlabel('时间 t');
ylabel('y(t)');
2. 刚性方程求解：ode15s
对于刚性方程，ode15s 更加合适。刚性方程指的是那些存在快速变化的解，会导致普通方法不稳定的问题。
示例：假设有以下刚性方程：
d
y
d
t
=
−
1000
⋅
(
y
−
cos
(
t
)
)
\frac{dy}{dt} = -1000 \cdot (y - \cos(t))
dtdy=−1000⋅(y−cos(t))
并且初始条件为 (
y
(
0
)
=
0
y(0) = 0
y(0)=0 )。
% 定义刚性微分方程
dydt = @(t, y) -1000 * (y - cos(t));
% 时间范围
tspan = [0, 0.1];
% 初始条件
y0 = 0;
% 使用 ode15s 求解刚性 ODE
[t, y] = ode15s(dydt, tspan, y0);
% 绘图
plot(t, y);
title('刚性方程解 \frac{dy}{dt} = -1000(y - \cos(t))');
xlabel('时间 t');
ylabel('y(t)');
3. 二阶常微分方程
二阶常微分方程可以通过将其转换为一组一阶方程来求解。例如，设有二阶方程：
d
2
y
d
t
2
+
5
d
y
d
t
+
6
y
=
0
\frac{d^2y}{dt^2} + 5\frac{dy}{dt} + 6y = 0
dt2d2y+5dtdy+6y=0
可以引入新的变量，将其转换为两个一阶方程：
- 设 ( y 1 = y y_1 = y y1=y )
- 设 ( y 2 = d y d t y_2 = \frac{dy}{dt} y2=dtdy )
则得到：
d
y
1
d
t
=
y
2
\frac{dy_1}{dt} = y_2
dtdy1=y2
d y 2 d t = − 6 y 1 − 5 y 2 \frac{dy_2}{dt} = -6y_1 - 5y_2 dtdy2=−6y1−5y2
MATLAB 代码如下：
% 定义微分方程组
dydt = @(t, y) [y(2); -6 * y(1) - 5 * y(2)];
% 时间范围
tspan = [0 5];
% 初始条件
y0 = [1; 0]; % y(0) = 1, y'(0) = 0
% 求解二阶常微分方程
[t, y] = ode45(dydt, tspan, y0);
% 绘图
plot(t, y(:,1)); % 只绘制 y1，即 y
title('二阶常微分方程解 \frac{d^2y}{dt^2} + 5\frac{dy}{dt} + 6y = 0');
xlabel('时间 t');
ylabel('y(t)');
4. 向量化求解
有时我们需要解多维度的常微分方程组，可以用向量化求解。例如下面的例子展示了一个两维系统的 ODE：
d
x
d
t
=
y
\frac{dx}{dt} = y
dtdx=y
d
y
d
t
=
−
x
\frac{dy}{dt} = -x
dtdy=−x
代码如下：
% 定义微分方程组
dydt = @(t, Y) [Y(2); -Y(1)];
% 时间范围
tspan = [0, 10];
% 初始条件
Y0 = [1; 0]; % x(0) = 1, y(0) = 0
% 求解 ODE
[t, Y] = ode45(dydt, tspan, Y0);
% 绘图
plot(t, Y(:,1), t, Y(:,2));
legend('x(t)', 'y(t)');
title('二元系统 \frac{dx}{dt} = y, \frac{dy}{dt} = -x');
xlabel('时间 t');
ylabel('解');
在 MATLAB 中，随机偏微分方程涉及的不仅是空间上的偏微分，还引入了时间和随机性。随机偏微分方程在许多领域都有应用，如金融、物理、气象学等，但其复杂性也更高。MATLAB 提供了一些工具来帮助处理和模拟简单的 SPDE，主要通过 有限差分、有限元等数值方法结合 随机过程 进行模拟。
下面介绍一些基本的随机偏微分方程模拟方法。
1. 随机扩散方程
扩散方程是最简单的偏微分方程之一，如果在其基础上引入噪声项，则可表示为：
∂
u
∂
t
=
D
∂
2
u
∂
x
2
+
σ
η
(
x
,
t
)
\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + \sigma \eta(x, t)
∂t∂u=D∂x2∂2u+ση(x,t)
其中，(
D
D
D ) 是扩散系数，(
σ
\sigma
σ ) 是噪声强度，(
η
(
x
,
t
)
\eta(x, t)
η(x,t) ) 是白噪声过程。
MATLAB 实现
% 设置参数
D = 0.1; % 扩散系数
sigma = 0.05; % 噪声强度
L = 1; % 空间长度
T = 1; % 时间长度
Nx = 100; % 空间网格数量
Nt = 500; % 时间步数
dx = L / (Nx - 1); % 空间步长
dt = T / Nt; % 时间步长
x = linspace(0, L, Nx); % 空间点
u = zeros(Nx, Nt); % 结果矩阵
% 初始条件
u(:, 1) = sin(pi * x); % 初始分布
% 模拟随机扩散过程
for n = 1:Nt-1
% 计算空间二阶导数
laplacian_u = ([u(2:end, n); 0] - 2 * u(:, n) + [0; u(1:end-1, n)]) / dx^2;
% 添加噪声项 (正态分布白噪声)
noise = sigma * sqrt(dt) * randn(Nx, 1);
% 更新 u 值
u(:, n+1) = u(:, n) + D * laplacian_u * dt + noise;
end
% 绘制结果
imagesc(linspace(0, T, Nt), x, u);
colorbar;
title('随机扩散方程解');
xlabel('时间 t');
ylabel('空间位置 x');
2. 二维空间中的扩散-反应方程
对于二维扩散-反应方程，可以引入随机项来模拟更复杂的过程，如模式形成。其一般形式为：
∂
u
∂
t
=
D
∇
2
u
+
f
(
u
,
v
)
+
σ
η
(
x
,
y
,
t
)
\frac{\partial u}{\partial t} = D \nabla^2 u + f(u, v) + \sigma \eta(x, y, t)
∂t∂u=D∇2u+f(u,v)+ση(x,y,t)
其中，(
f
(
u
,
v
)
f(u, v)
f(u,v) ) 是反应函数，通常代表某种动力学过程。
MATLAB 实现（简单随机项）
% 设置参数
D = 0.1; % 扩散系数
sigma = 0.01; % 噪声强度
Lx = 1; % x 方向长度
Ly = 1; % y 方向长度
T = 1; % 总时间
Nx = 50; % x 方向网格数
Ny = 50; % y 方向网格数
Nt = 100; % 时间步数
dx = Lx / (Nx - 1); % x 方向步长
dy = Ly / (Ny - 1); % y 方向步长
dt = T / Nt; % 时间步长
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
% 初始条件
u = zeros(Nx, Ny, Nt);
u(:, :, 1) = sin(pi * x') * sin(pi * y);
% 时间演化
for n = 1:Nt-1
% 计算二维拉普拉斯项
laplacian_u = (circshift(u(:, :, n), [1, 0]) + circshift(u(:, :, n), [-1, 0]) ...
+ circshift(u(:, :, n), [0, 1]) + circshift(u(:, :, n), [0, -1]) ...
- 4 * u(:, :, n)) / (dx^2);
% 随机噪声项
noise = sigma * sqrt(dt) * randn(Nx, Ny);
% 更新 u 值
u(:, :, n+1) = u(:, :, n) + D * laplacian_u * dt + noise;
end
% 可视化随机扩散-反应结果
for n = 1:Nt
imagesc(x, y, u(:, :, n));
colorbar;
title(['时间 t = ' num2str(n * dt)]);
xlabel('x');
ylabel('y');
pause(0.1);
end
3. 其他随机偏微分方程求解方法
MATLAB 直接求解随机偏微分方程时，有时需要使用有限元方法、随机数值积分等更高级的技术。为了简化，我们可以基于分块处理的思路，将空间和时间离散化，再引入随机项来求解。MATLAB 的 PDE 工具箱也可以在扩展后用于随机偏微分方程，尽管 MATLAB 的基础函数并不直接支持随机偏微分方程。
例如，可以通过 MATLAB 的 PDE 工具箱配合自定义噪声项来模拟简单的随机扩散过程。以上代码示例展示了 MATLAB 中随机偏微分方程的基本模拟方法。虽然 MATLAB 没有专用的随机偏微分方程求解工具，但可以通过引入随机噪声项并结合差分或有限元方法来模拟。
