How to plot Ramsey Phase Plane in Matlab

7 次查看(过去 30 天)
Hey,
I am currently working on a script meant to plot the trajectory, saddle paths and the phase plane of the basic RKC model. However I continue to get the following error: Error using streamslice. U,V must all be size 2x2 or greater. Error in Code_Diss (line 96). h = streamslice(K, C, dK, dC, 2, 'noarrows', 'cubic');
I havent worked much in Matlab yet so forgive me if this is a very obvious error.
Thanks!
clear; % clear Workspace
clc; % clear Command Window
close all; % close Figures
disp('----------------------------------');
disp(' Ramsey-Cass-Koopmans model: ');
disp('----------------------------------');
%% [0] define variables
global alpha delta rho n g theta kss css k0
syms t L(t) K(t) C(t) a(t) r(t) w(t) u(t) rho n theta
L(t) = exp(n.*t); % total population at t
n = diff(L(t),t)/L(t); % labor growth rate
c(t) = C(t)/L(t); % per capita consumption
k(t) = K(t)/L(t); % per capita capital
%% [1] define parameters
rho = 0.04; % discount rate
theta = 1.5; % IES
nu = 0.35;
sigma = 1;
delta = 0.05; % depreciation rate
n = 0.01; % population growth rate
alpha = 0.5 ; % capital share
c_0 = 3 ; % Initial per capital consumption
k_0 = 5 ; % Initial per capita capital (wealth)
initial = [3, 5]; % c(0) and k(0)
kmin=0; % smallest capital stock to consider
kmax=100; % largest capital stock to consider
u = ((C(t)/L(t))^(1-theta))/(1-theta) ; % utility at time t
simplify(int(u*exp(-(rho-n)*t),0,inf)); % lifetime utility
%% [2] compute steady state for c and k
syms k l c ;
f = k.^ (alpha) * l.^(1 - alpha) ;
y = f
df_dk = alpha * k.^(alpha - 1) * l.^(1 - alpha) ;
df_dl = (1 - alpha) * k.^(alpha) * l^(- alpha) ;
% Write down the equations for dk/dt and dc/dt.
k_dot = y - c - (delta - n) * k == 0 ;
c_dot = c * (((df_dk - delta - rho)/ theta )+ n) == 0 ;
% Steady state value of k obtained when dc/dt = 0.
kss = df_dk - delta - rho - ( n * theta ) ;
% Compute the steady-state curve c = c*(k) obtained when dk/dt = 0.
css = y - (delta + n) * k ;
%% Define sample values for the x-axis (per-capita wealth).
k_root = (n + delta) ^ (1 / (alpha - 1));
nPoints = 50;
k_init = 0;
k_fine = linspace(k_init, k_root, nPoints);
%% Define sample values for the y-axis (per-capita consumption).
c_fine = linspace(k_init, c_0, nPoints);
%% Create a lattice of (k, c) points for the stream slices.
[K, C] = meshgrid(k_fine, c_fine);
dK = k_dot
dC = c_dot
kss = 4
css = 3
disp('steady-state:')
fprintf('k* = %14.6f \n',kss);
fprintf('c* = %14.6f \n',css);
%plot_ramsey
%% Visualize the phase portrait.
figure('Units', 'Normalized', 'Position', 0.25*[1, 1, 2, 2])
h = streamslice(K, C, dK, dC, 2, 'noarrows', 'cubic');
set(h, 'LineStyle', ':', 'LineWidth', 1.5)
% Hide all but one streamslice from the legend.
set(h(2:end), 'HandleVisibility', 'off')
% Plot the steady-state vertical line and smooth curve.
hold on
plot(k_fine, c_star, 'r', 'LineWidth', 2)
plot([k_steady, k_steady], ylim(), 'm', 'LineWidth', 2)
xlabel('{\it{k}} (capital per capita)')
ylabel('{\it{c}} (consumption per capita)')
title('Ramsey-Cass-Koopmans Capital-Consumption Phase Plane')
grid on
% Increase the font size.
ax = gca;
ax.FontSize = 15;
% Add the legend.
legend({'Stream slices', '$dk/dt = 0$', '$dc/dt = 0$'}, ...
'Interpreter', 'LaTeX', ...
'FontSize', 15)
pause
close

回答(1 个)

SAI SRUJAN
SAI SRUJAN 2023-10-5
Hi Sabrina Haumann,
I understand that you are attempting to plot the trajectory, saddle paths, and phase plane, for the fundamental RKC model. However, you are encountering an issue when utilizing the streamslice function in MATLAB.
In the line 96,
%streamslice(X,Y,U,V);
h=streamslice(K,C,dK,dC,2,"noarrows","cubic");
The aforementioned code line encounters an error due to the requirement of "dK" and "dC" in streamslice to have a size of 2x2 or greater. The acceptable components, "dK" and "dC", should be either a 2D array or a 3D array.
However, in your code, "dK" and "dC" have a size of 1x1 syms. Consequently, the function "streamslice" fails to execute successfully.
You may find valuable insights in the provided documentation, which includes examples demonstrating the proper usage of streamslice.

类别

Help CenterFile Exchange 中查找有关 Vector Fields 的更多信息

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by