基于Matlab的离散化SDE系统多轨迹随机模拟技术需求

Got it, let's walk through implementing this SDE simulation in MATLAB step by step. I'll break this into clear, actionable chunks so you can replicate it easily.

1. Define Core Parameters #

First, let's lock in all the fixed values we need upfront—adjust these to match your problem's specifics:

• Number of trajectories: num_trajectories = 10000
• Total simulation days: num_days = 100
• Daily time step (since we're simulating day-by-day, this is 1): dt = 1
• Constant c (fill in your actual value here): c = 1; % Example placeholder
• Initial conditions for x and y:

  x0 = 0.5; % Example starting value for x
  y0 = 1.0; % Example starting value for y
  

2. Weekday/Weekend Logic #

We need to switch f1 and f2 behavior between weekdays (Mon-Fri) and weekends (Sat-Sun). A simple, date-independent way to do this uses modular arithmetic:

% Precompute a logical array: true = weekday, false = weekend
is_weekday = arrayfun(@(t) mod(t-1,7) < 5, 1:num_days);

This works because every 7-day cycle, days 1-5 count as weekdays, and days 6-7 as weekends.

3. Initialize Simulation Arrays #

Let's set up arrays to store the full trajectory of x and y for all 10000 paths:

% Rows = trajectories, columns = time steps (days 0 to 100)
x = zeros(num_trajectories, num_days + 1);
y = zeros(num_trajectories, num_days + 1);

% Set initial conditions for all trajectories
x(:,1) = x0;
y(:,1) = y0;

4. Run the Stochastic Simulation #

Now we'll loop through each day, compute f1/f2 based on the day type, and update x/y using your SDE discretization. We'll add a tiny epsilon to avoid division by zero (critical since x and y appear in denominators):

% Small value to prevent division by zero
eps_val = 1e-8;

for t = 1:num_days
    % Grab current x and y values for all trajectories
    x_t = x(:,t);
    y_t = y(:,t);
    
    % Ensure no zero values in denominators
    x_t(x_t < eps_val) = eps_val;
    y_t(y_t < eps_val) = eps_val;
    
    % Calculate f1 and f2 based on weekday/weekend
    if is_weekday(t)
        % Monday-Friday rules
        f1 = 2 .* x_t.^2 ./ (y_t + x_t + c);
        f2 = y_t.^2;
    else
        % Saturday-Sunday rules (note: your original question cut off f2 for weekends—replace the placeholder below!)
        f1 = x_t ./ y_t;
        f2 = % Insert your weekend f2 formula here
    end
    
    % Generate standard normal random numbers for each trajectory
    rand1 = randn(num_trajectories, 1);
    rand2 = randn(num_trajectories, 1);
    
    % Update x and y using the given SDE equations
    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

Important: Don't forget to fill in the missing weekend f2 formula from your original problem statement!

5. Optional: Visualize Results #

Plotting a subset of trajectories (or the mean across all paths) helps verify the simulation works:

% Plot first 10 trajectories of x
figure;
plot(0:num_days, x(1:10,:)');
xlabel('Day');
ylabel('x(t)');
title('Sample Trajectories of x(t)');
legend('Traj 1','Traj 2','Traj 3','Traj 4','Traj 5',...
       'Traj 6','Traj 7','Traj 8','Traj 9','Traj 10');
grid on;

% Plot mean x value across all 10000 trajectories
figure;
plot(0:num_days, mean(x,1));
xlabel('Day');
ylabel('Mean x(t)');
title('Average Trajectory of x(t)');
grid on;

Quick Tips #

• Reproducibility: Add rng(42) (or any seed number) at the start to fix the random number generator and get the same results every run.
• Stability: If you notice extreme values, tweak the initial conditions or eps_val to keep x/y in a reasonable range.

内容的提问来源于stack exchange，提问作者maryam john