I am trying to solve the system of two linear differential equations and create a phase diagram to asses the stability of the system.

\dot{\dot{y}\ =-\delta\gamma\beta(y-y_n)-\delta\gamma\lambda(p-p_t)}

\dot{\dot{p}=\alpha(y-y_n)}

Answer 1

[Sam Chak]

Hi @Chris Herdelin

If you are referring to the phase diagram as the 2-D phase plane plot, you can use a special output function called 'OutputFcn' and specify the function handle as '@odephas2' to plot it.

By the way, there appear to be errors (extra \dot) in the LaTeX typesetting. The example below uses this system:

Please copy and paste the provided code into a script and save it as "mySystem.m" in the current working folder. To open the script, type "edit mySystem" in the Command Window, and to run it, simply enter "mySystem" in the Command Window.

tspan   = [0 10];
p0      =  1;
y0      =  0;
x0      = [p0; y0];
options = odeset('OutputFcn', 'odephas2');
[t, x]  = ode45(@odefcn, tspan, x0, options);

fig                 = gcf;
ax                  = fig.CurrentAxes;
ax.XGrid            = "On";
ax.YGrid            = "On";
ax.XLim             = [-1.0, 1.0];
ax.YLim             = [-0.5, 0.5];
ax.XLabel.String    = 'p(t)';
ax.YLabel.String    = 'y(t)';
ax.Title.String     = 'Phase Plane Plot';

%% System of two linear differential equations
function dx = odefcn(t, x)
    % definitions
    p       = x(1);
    y       = x(2);

    % parameters
    alpha   = 1;
    beta    = 2;
    gamma   = 1;
    delta   = 1;
    lambda  = 1;
    pt      = 0;
    yn      = 0;
    
    % linear ODEs
    dx(1,1) = alpha*y;
    dx(2,1) = - delta*gamma*beta*(y - yn) - delta*gamma*lambda*(p - pt);
end

[Chris Herdelin]

Sam this looks great! How do I get it to continue to spiral more so that I can show that the system is stable i.e. tending toward zero?

[Chris Herdelin]

Does it matter that I have Matlab 2023A? Should I update to the most revent version? Thanks again!

[Chris Herdelin]

BTW- You are correct regarding the extra dot in the LaTex code...

[Sam Chak]

Hi @Chris Herdelin

There is no need to upgrade the version. If you wish to observe a more pronounced effect of the spiral trajectory, you will need to adjust the parameters. In particular, the value of lambda should be significantly greater than beta.

Please let me know if the examples are satisfactory to you.

tspan   = [0 10];
p0      =  1;   % initial value of state variable x1
y0      =  0;   % initial value of state variable x2
x0      = [p0; y0];
[t, x]  = ode45(@odefcn, tspan, x0);

%% Plot x2 (y-axis) vs. x1 (x-axis)
plot(x(:,1), x(:,2)), grid on, 
axis([-1 1 -1 1])
xlabel('p(t)'), 
ylabel('y(t)'), 
title('Phase Plane Plot')

%% System of two linear differential equations
function dx = odefcn(t, x)
    % definitions
    p       = x(1);
    y       = x(2);

    % parameters
    alpha   = 1;
    beta    = 1;
    gamma   = 1;
    delta   = 1;
    lambda  = 2;
    pt      = 0;
    yn      = 0;
    
    % linear ODEs
    dx(1,1) = alpha*y;
    dx(2,1) = - delta*gamma*beta*(y - yn) - delta*gamma*lambda*(p - pt);
end

[Sam Chak]

The syntax is "plot(x_axis, y_axis)". First argument is x-axis signal, and second argument is y-axis signal.

Answer 2

[Joshua Levin Kurniawan]

First, you can define the ODE function

function dydt = myODE(t, y, gamma, beta, yn, delta, lambda, pt, alpha)
    y1 = y(1); 
    y2 = y(2);
    
    dy1dt = y2;
    dy2dt = -delta*gamma*beta*(y1-yn) - delta*gamma*lambda*(y2-pt);
    dydt = [dy1dt; dy2dt];
end

For example:

gamma = 1;
beta = 1;
yn = 0;
delta = 1;
lambda = 1;
pt = 0;
alpha = 1;

% Initial conditions
y0 = [0; 0]; % [y; p] initial condition

% Simulation time soan (in seconds)
tspan = [0 10];

% Solve the differential equations
[t, y] = ode45(@(t, y) myODE(t, y, gamma, beta, yn, delta, lambda, pt, alpha), tspan, y0);

% Plotting result
plot(t, y(:, 1), 'b', t, y(:, 2), 'r');
xlabel('Time');
ylabel('y and p');
legend('y', 'p');

If you want to see the Bode diagram of the system, you can transform it to Laplace domain transfer function:

num_y = [-delta*gamma*beta*yn - delta*gamma*lambda*pt];
den_y = [1, 0, delta*gamma*beta];
G_y = tf(num_y, den_y);

num_p = [-alpha*yn];
den_p = [1, 0, -alpha];
G_p = tf(num_p, den_p);

% Create bode diagram
bode(G_y);
bode(G_p);

[Chris Herdelin]

Haven't tried this yet, but wow this looks promising! Thanks so much! I will let you know if I hae any issues running the code!

[Chris Herdelin]

[Torsten]

Instead of "untitled", save your function as "myOde.m" in your working directory.

[Chris Herdelin]

Any change you show how you are running the script?

[Torsten]

Which MATLAB version do you use ?

Under the recent version, it works (at least technically).

gamma = 1;
beta = 1;
yn = 0;
delta = 1;
lambda = 1;
pt = 0;
alpha = 1;
% Initial conditions
y0 = [0; 0]; % [y; p] initial condition
% Simulation time soan (in seconds)
tspan = [0 10];
% Solve the differential equations
[t, y] = ode45(@(t, y) myODE(t, y, gamma, beta, yn, delta, lambda, pt, alpha), tspan, y0);
% Plotting result
plot(t, y(:, 1), 'b', t, y(:, 2), 'r');
xlabel('Time');
ylabel('y and p');
legend('y', 'p');

num_y = [-delta*gamma*beta*yn - delta*gamma*lambda*pt];
den_y = [1, 0, delta*gamma*beta];
G_y = tf(num_y, den_y);
num_p = [-alpha*yn];
den_p = [1, 0, -alpha];
G_p = tf(num_p, den_p);
% Create bode diagram
bode(G_y);

bode(G_p);

function dydt = myODE(t, y, gamma, beta, yn, delta, lambda, pt, alpha)
    y1 = y(1); 
    y2 = y(2);
    
    dy1dt = y2;
    dy2dt = -delta*gamma*beta*(y1-yn) - delta*gamma*lambda*(y2-pt);
    dydt = [dy1dt; dy2dt];
end

[Sam Chak]

Hi @Chris Herdelin,

It is not advisable to manually enter multiple parameters one by one and do 'long' function calls in the Command Window. This is because you would get into a hassle of re-entering them each time you want to run the simulation, especially if the variables in the Workspace are cleared (either by you or after exiting MATLAB).

If you do not wish to create and run a MATLAB script (.m or .mlx), it would be best to store the parameters in a simple Notepad file on your Desktop. This way, you can easily locate, access and copy/paste the parameters and commands whenever you need them.