ODE simulation: issue with the size of arrays

2 次查看（过去 30 天）
Alessandro Bellocchi
编辑： Davide Masiello 2023-4-13
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');
  1 个评论
Davide Masiello
Davide Masiello 2023-4-13
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Davide Masiello
Davide Masiello 2023-4-13
编辑：Davide Masiello 2023-4-13
Ran in:
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')

类别

在 Help Center 和 File Exchange 中查找有关 Interest-Rate Instruments 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by