Matlab 模拟（随机）答案

【问题标题】：Matlab Simulation(Stochastic)Matlab 模拟（随机）
【发布时间】：2018-05-07 12:32:46
【问题描述】：

假设我有一个形式为 SDE 的离散系统

x(:, t+1) = x(:, t) + f1(x(:, t)).*x(:, t)*dt + f2(x(:, t))./(x(:, t).*y(:, t))* sqrt(dt)*rand1;

y(:, t+1) = f2(x(:, t)).*y(:, t)./x(:, t)*dt + f1(x(:, t)).*y(:, t)*sqrt(dt)*rand2;

我想用 10000 条轨迹模拟系统，

对于时间 t = 100 天，例如：从星期一到星期五，

f1(x(:, t)) = 2*x(:, t).^2./(y(:, t) + x(:, t) + c)，和

f2(x(:, t)) = y(:, t).^2; 而周六和周日

f1(x(:, t)) = x(:, t)./y(:, t) 和 f2(x(:, t)) = y(:, t); 如何模拟 SDE 系统？

这是我的方法

dt = 0.01;
time = 100;
num_iteration = ceil(time / dt);
num_trajectory = 10000; 
%% Initial Values
y0 = 1; 
x0 = 1;
y = zeros(num_trajectory, num_iteration) + y0; 
x = zeros(num_trajectory, num_iteration) + x0; 
days = 0;

for t=1: num_iteration
    current_time = t * dt;
    rand1 = randn(num_trajectory, 1);
    rand2 = randn(num_trajectory, 1);

    if ceil(current_time) == current_time
        days = days+1;

        if (mod(days, 7) | mod(days+1, 7)) == 0
            f1 = 2*x(:, t).^2./(y(:, t) + x(:, t) + c);
            f2 = y(:, t).^2;
        else
            f1 = x(:, t)./y(:, t);
            f2 = y(:, t); 
        end
    end

    x(:, t+1) = x(:, t) + f1*x(:, t)*dt + f2/(x(:, t).*y(:, t))* sqrt(dt)*rand1;
    y(:, t+1) = f2*y(:, t)./x(:, t)*dt + f1*y(:, t)*sqrt(dt)*rand2;   
end

【问题讨论】：

文本不包含实际问题。你想达到什么目标，你已经做了什么？
对于给定的dt，并假设变量中的行代表 10000 个“轨迹”，那么对于每一行，您需要一个合适的代表 t=1 的“起点”，您将使用 x(:,t) 来查找 x(:,t+1)，即 x 表示 t=2。让 t=2 你继续找到 t=3 等的 x，直到 t = 100。rand 函数涉及的事实，在每一步引入一些随机“噪音”，这意味着即使起点所有行都相同（例如，所有行的 x(t=0) = 0），每行遵循的最终轨迹将不同。
到目前为止，这是我所做的：
@TasosPapastylianou，实际上系统比我发布的要复杂得多。我做了你提到的一切，但我仍然没有得到预期的结果。

标签： matlab simulation stochastic

【解决方案1】：

你的方法看起来不错。不过，您的代码中有一个逻辑错误。在行中

if (mod(days, 7) | mod(days+1, 7)) == 0

表达式 (mod(days, 7) | mod(days+1, 7)) 将始终计算为 1（尝试找出原因），因此 (mod(days, 7) | mod(days+1, 7)) == 0 将始终为 false，并且您的 if 语句将始终将控制权传递给 else 部分。

因此应该改为

if mod(days, 7) == 0 || mod(days+1, 7) == 0

但这也令人困惑（而且您还没有在代码中记录“0”是星期几）。

更清楚的是：

if (
  mod (days, 7) == 0 % day is a sunday
  || 
  mod (days, 7) == 6 % day is a saturday
)
    % do stuff
else
    % do other stuff
end

更好的是，创建一个小函数 isWeekend 为您执行该测试，从而生成超级清晰的代码，例如

if isWeekend(days)
  % do stuff
else 
  % do other stuff
end

【讨论】：

非常感谢，你们的cmets很有帮助。