Matlab 仿真(随机)
Matlab Simulation(Stochastic)
假设我有一个形式为
的 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
你的方法看起来不错。但是,您的代码中存在逻辑错误。在行
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
假设我有一个形式为
的 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
你的方法看起来不错。但是,您的代码中存在逻辑错误。在行
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