Matlab Simulation (Стохастик)

Question

Предполагая, что у меня есть дискретизированная система 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

Answer 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 всегда будет ложным, а оператор 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