ODE simulation: issue with the size of arrays

Hello, everyone. I am trying to simulate the evolution of the following model. However, I get an error with the size of the arrays (Unable to perform assignment because the size of the left side is 1-by-2 and the size of the right side is 1-by-330201). Can anyone help me solve it?

% Define the parameters
m = 0.25;
beta = 1.25;
rho = 0.04;
delta = 0.03;
alpha = 0.2;
pphi = 0.5;
mu = 0.2;
gamma = 1;

% Define the initial condition
phi0 = 10;
theta0 = 1;

% Define the values of sigma for the two situations
sigma_values = [0.1, 0.5];

% Define the time span for the simulation
tspan = [0 500];

% Define the Wiener process
dW = makedist('Normal', 'mu', 0, 'sigma', sqrt(diff(tspan)));

% Define the number of simulations
num_simulations = 5;

% Initialize arrays to store the phi values for each simulation
phi_values = zeros(length(sigma_values), length(tspan), num_simulations);
average_phi_values = zeros(length(sigma_values), length(tspan));

% Perform multiple simulations and store the phi values
for sim = 1:num_simulations
    for i = 1:length(sigma_values)
        % Define the current value of sigma
        sigma = sigma_values(i);
        
        % Define the function to be solved
        finance = @(t, y) [
            (1/gamma)*((m*y(2)^beta)/(rho+delta-(alpha*(beta*sigma)^2)/2)-(pphi-mu)-delta*y(1));
            sigma*y(2)*dW.random()];
        
        % Solve the differential equation numerically
        [t, y] = ode45(finance, tspan, [phi0, theta0]);
        
        % Store the phi values for the current simulation and sigma value
        phi_values(i, :, sim) = y(:, 1)';
    end
end

% Calculate the average phi value over all simulations for each time point and sigma value
for i = 1:length(sigma_values)
    for j = 1:length(tspan)
        average_phi_values(i, j) = mean(phi_values(i, j, :));
    end
end

% Plot the average phi values against time for each sigma value
hold on
for i = 1:length(sigma_values)
    plot(tspan, average_phi_values(i, :), 'LineWidth', 2);
end

% Add legend, labels, and title
legend(['\sigma = ', num2str(sigma_values(1))], ['\sigma = ', num2str(sigma_values(2))]);
xlabel('Time');
ylabel('Average \phi_t');
title('Evolution of Finance over Time for Different Values of \sigma and Wiener Process');

[Davide Masiello]

Unfortunately, you do not know the size of the t vector before the ODE integration itself.

Hence, you cannot initialize phi_values the way you did.

Even if you didn't initialize it, you'd probably end up with a size error, because each simulation will require a different number of steps to complete the integration.

This is due to the fact that MatLab ode solvers are adaptive solvers, with variable step size.

Answer 1

[Davide Masiello]

In addition to my comment above, here's a solution to your problem.

It's based on the fact that you can specify tspan as a vector of points at which the solution must be reported.

This provides a homogeneous output regardless of the choice of the integration time steps.

clear,clc
% Define the parameters
m = 0.25;
beta = 1.25;
rho = 0.04;
delta = 0.03;
alpha = 0.2;
pphi = 0.5;
mu = 0.2;
gamma = 1;
% Define the initial condition
phi0 = 10;
theta0 = 1;
% Define the values of sigma for the two situations
sigma_values = [0.1, 0.5];
% Define the time span for the simulation
tspan = linspace(0,500,1000);
% Define the Wiener process
dW = makedist('Normal', 'mu', 0, 'sigma', sqrt(tspan(end)-tspan(1)));
% Define the number of simulations
num_simulations = 5;
% Initialize arrays to store the phi values for each simulation
phi_values = zeros(length(sigma_values), length(tspan), num_simulations);
average_phi_values = zeros(length(sigma_values), length(tspan));
% Perform multiple simulations and store the phi values
for sim = 1:num_simulations
    for i = 1:length(sigma_values)
        % Define the current value of sigma
        sigma = sigma_values(i);
        % Define the function to be solved
        finance = @(t, y) [
            (1/gamma)*((m*y(2)^beta)/(rho+delta-(alpha*(beta*sigma)^2)/2)-(pphi-mu)-delta*y(1));
            sigma*y(2)*dW.random()];
        % Solve the differential equation numerically
        [t, y] = ode45(finance, tspan, [phi0, theta0]);
        % Store the phi values for the current simulation and sigma value
        phi_values(i, :, sim) = y(:, 1)';
    end
end
% Calculate the average phi value over all simulations for each time point and sigma value
for i = 1:length(sigma_values)
    for j = 1:length(tspan)
        average_phi_values(i, j) = mean(phi_values(i, j, :));
    end
end
% Plot the average phi values against time for each sigma value
plot(tspan, average_phi_values, 'LineWidth', 2);

Warning: Imaginary parts of complex X and/or Y arguments ignored.

% Add legend, labels, and title
legend(['\sigma = ', num2str(sigma_values(1))], ['\sigma = ', num2str(sigma_values(2))]);
xlabel('Time');
ylabel('Average \phi_t');
title('Evolution of Finance over Time for Different Values of \sigma and Wiener Process')