Does anyone know how to use an IF statement with ODE's? I am trying to create an IF statement sigma > 0 , tau > 0 if dx(3,1) > 0, dx(1,1) > 0 and sigma = 0 , tau = 0 otherwise

close all
clear all
clc
tspan   = [0 100];
y0      =  0.7;   % initial value of state variable x1
r0      =  0.7;   % initial value of state variable x2
w0      =  0.8;   % initial value of state variable x3
x0      = [y0; r0; w0];
[t, x]  = ode45(@odefcn, tspan, x0);
%% Plot x2 (y-axis) vs. x1 (x-axis)
subplot(3,1,1);
plot(t,x(:,1));
xlabel('Time')
ylabel('Output');
subplot(3,1,2);
plot(t,x(:,2));
xlabel('Time')
ylabel('Policy Rate');
subplot(3,1,3);
plot(t,x(:,3));
xlabel('Time')
ylabel('Wage Share');
y0 = [0.07; 0.07; 0.8];
%% System of three linear differential equations
function dx = odefcn(t, x)
% definitions
y       = x(1);
r       = x(2);
w       = x(3);

% parameters
alpha   = 1.0;
beta    = 1.0;
gamma   = 0.3;
delta   = 0.3;
mu      = 0.4;
lambda  = 1.0;
theta   = 0.4;
omega   = 0.4;
sigma   = 0.1;
tau     = 0.1;

% linear ODEs
dx(1,1) = alpha * y - beta * r * y;
dx(2,1) = - omega * r + gamma * w + sigma * w + delta * y + tau * y;
dx(3,1) = - theta * w + lambda * y * w - mu * w * w;
end

Answer 1

[Sam Chak]

Hi @Christopher

I simulated 5 cases below, including the interpretations of @Torsten and @William Rose, while keeping the settings for both relative error and absolute error tolerances at their default values.

• Case 1: and
• Case 2: and
• Case 3: σ and τ continuously cycle between 0 and 1 (2nd derivative exists)
• Case 4: Torsten's interpretation
• Case 5: William's interpretation (2nd derivative exists)

According to the theoretical analysis, Cases 4 and 5 should exhibit oscillatory behavior, similar to Cases 1, 2, and 3. This is due to the coefficients of the states w and y in the state equation dx(2) remaining significantly positive, regardless of the abrupt changes in sigma and tau.

tspan   = [0 100];
y0      =  0.7;     % initial value of state variable x1
r0      =  0.7;     % initial value of state variable x2
w0      =  0.8;     % initial value of state variable x3
x0      = [y0; r0; w0];
% opt     = odeset('Reltol', 1e-6, 'Abstol', 1e-12);

cs      = 1:5;      % Choose Case from 1 to 5

for i = 1:numel(cs)
    [t, x]  = ode45(@(t, x) odefcn(t, x, cs(i)), tspan, x0);

    dx1 = zeros(1, numel(t));
    dx3 = zeros(1, numel(t));
    for j = 1:numel(t)
        [~, dx1(j), dx3(j)] = odefcn(t(j), x(j,:).', cs(i));
    end
    
    %% Plot results
    figure
    subplot(3,1,1);
    plot(t, x(:,1)); grid on
    xlabel('Time'), ylabel('Output');
    
    subplot(3,1,2);
    plot(t, x(:,2)); grid on
    xlabel('Time'), ylabel('Policy Rate');
    
    subplot(3,1,3);
    plot(t, x(:,3)); grid on
    xlabel('Time'), ylabel('Wage Share');
    
    sgtitle(sprintf('Case %d', cs(i)))
    
    figure
    subplot(211)
    plot(t, dx1), grid on
    xlabel('Time'), ylabel('dx1'), title('Time response of dx1')
    
    subplot(212)
    plot(t, dx3), grid on
    xlabel('Time'), ylabel('dx3'), title('Time response of dx3')
    
    sgtitle(sprintf('Case %d', cs(i)))
end

●

%% System of three linear differential equations
function [dx, dx1, dx3] = odefcn(t, x, cs)
    % definitions
    y       = x(1);
    r       = x(2);
    w       = x(3);

    % parameters
    alpha   = 1.0;
    beta    = 1.0;
    gamma   = 0.3;
    delta   = 0.3;
    mu      = 0.4;
    lambda  = 1.0;
    theta   = 0.4;
    omega   = 0.4;
    sigma   = 0.1;
    tau     = 0.1;
    
    dx1     =   alpha*y -   beta*r*y;
    dx3     = - theta*w + lambda*y*w - mu*w*w;

    % ODEs
    dx(1,1) =   alpha*y -   beta*r*y;
    
    switch cs
        case 1
            %% Situation 1: sigma ~= 0 && tau ~= 0
            dx(2,1) = - omega*r + (gamma + sigma)*w + (delta + tau)*y;
        case 2
            %% Situation 2: sigma  = 0 && tau  = 0
            dx(2,1) = - omega*r + gamma*w + delta*y;
        case 3
            %% Situation 3: sigma & tau cycle between 0 and 1 (2nd derivative exists)
            Tp      = 1e-2;     % Time period in the one cycle
            cycle   = (1 + sin(2*pi/Tp*t))/2;
            dx(2,1) = - omega*r +  (gamma + sigma*cycle)*w + (delta + tau*cycle)*y; 
        case 4
            %% Situation 4: Sam's reinterpretation of Torsten's interpretation
            dx(2,1) = - omega*r +  (gamma + sigma*(dx1 > 0 && dx3 > 0))*w + (delta + tau*(dx1 > 0 && dx3 > 0))*y;
        case 5
            %% Situation 5: William's interpretation
            ws      = 0.01;     % transition width for sigma (units: same as sigma)
            wt      = 0.01;     % transition width for tau (units: same as tau)
            dx(2,1) = - omega*r +  (gamma + sigma/(1 + exp(-dx3/ws)))*w + (delta + tau/(1 + exp(-dx1/wt)))*y;    
        otherwise
            warning('!!!')
    end

    dx(3,1) = - theta*w + lambda*y*w - mu*w*w;
end

[Christopher]

@Sam Chak This is amazing! The oscilations are exactly what I would hope to see! I am going to work through each case and give an update once I settle on which case makes the most sense! It looksl like case three has a lot of potential, but I have to look at a phase plane plot to make the ultimtate decision.

[Christopher]

@Sam Chak Thanks again for the updated script! I was able to run it fully, but how do I just isolate case three?

[Walter Roberson]

To isolate case three use

cs = 3;

instead of

cs = 1:5;

At some point you would remove all the "case" except for 3 and "otherwise". After that you would remove cs as an input and remove the "switch" and "case 3" and "otherwise".

[Christopher]

@Walter Roberson That is very helpful!

[Christopher]

@Sam Chak I was hoping to see a declining trend in the wage share and output when looking at the 2d Phase diagram, but it doesn't look like that is the result. Could I be missing something?

[Sam Chak]

Hi @Christopher

If you run the simulation for a longer duration, you will observe that the responses of Case 3 behave like chaotic systems. In fact, the phase portrait resembles Jupiter's Great Red Spot. You need to introduce an 'actuation' variable that can steer the responses toward your desired outcome.

However, at this moment, I can only identify the coefficients as the means to influence the responses. This is analogous to adjusting the damping ratio and stiffness in an unforced mass-spring-damper system.

% create color map (cmap)
cols        = 2160;
rows        = 3840;
top_color   = [255,  89,  89]/255;  % Sunset Orange (top color)
bot_color   = [255, 209, 117]/255;  % Topaz (bottom color)
cmap        = interp1([0, 1], [bot_color; top_color], linspace(0, 1, cols));

tspan   = 0:0.1:10000;
y0      =  0.7;     % initial value of state variable x1
r0      =  0.7;     % initial value of state variable x2
w0      =  0.8;     % initial value of state variable x3
x0      = [y0; r0; w0];
opt     = odeset('Reltol', 1e-6, 'Abstol', 1e-12);

cs      = 3;      % Choose Case from 1 to 5

for i = 1:numel(cs)
    [t, x]  = ode45(@(t, x) odefcn(t, x, cs(i)), tspan, x0);

    dx1 = zeros(1, numel(t));
    dx3 = zeros(1, numel(t));
    for j = 1:numel(t)
        [~, dx1(j), dx3(j)] = odefcn(t(j), x(j,:).', cs(i));
    end
    
    %% Plot results
    figure
    subplot(3,1,1);
    plot(t, x(:,1)); grid on
    xlabel('Time'), ylabel('Output');
    
    subplot(3,1,2);
    plot(t, x(:,2)); grid on
    xlabel('Time'), ylabel('Policy Rate');
    
    subplot(3,1,3);
    plot(t, x(:,3)); grid on
    xlabel('Time'), ylabel('Wage Share');
    
    sgtitle(sprintf('Case %d', cs(i)))
    
    figure
    % plot(x(:,1), x(:,3)), grid on
    z   = x(:,3)';
    patch([x(:,1)' NaN], [x(:,3)' NaN], [z NaN], [z NaN], 'EdgeColor', 'interp'); 
    colormap(cmap);
    grid on
    xlabel('Output'), ylabel('Wage Share');
    title('Phase Portrait of Wage Share vs. Output')
end

%% System of three linear differential equations
function [dx, dx1, dx3] = odefcn(t, x, cs)
    % definitions
    y       = x(1);
    r       = x(2);
    w       = x(3);

    % parameters
    alpha   = 1.0;
    beta    = 1.0;
    gamma   = 0.3;
    delta   = 0.3;
    mu      = 0.4;
    lambda  = 1.0;
    theta   = 0.4;
    omega   = 0.4;
    sigma   = 0.1;
    tau     = 0.1;
    
    dx1     =   alpha*y -   beta*r*y;
    dx3     = - theta*w + lambda*y*w - mu*w*w;

    % ODEs
    dx(1,1) =   alpha*y -   beta*r*y;
    
    switch cs
        case 1
            %% Situation 1: sigma ~= 0 && tau ~= 0
            dx(2,1) = - omega*r + (gamma + sigma)*w + (delta + tau)*y;
        case 2
            %% Situation 2: sigma  = 0 && tau  = 0
            dx(2,1) = - omega*r + gamma*w + delta*y;
        case 3
            %% Situation 3: sigma & tau cycle between 0 and 1 (2nd derivative exists)
            Tp      = 1e-2;     % Time period in the one cycle
            cycle   = (1 + sin(2*pi/Tp*t))/2;
            dx(2,1) = - omega*r +  (gamma + sigma*cycle)*w + (delta + tau*cycle)*y; 
        case 4
            %% Situation 4: Sam's reinterpretation of Torsten's interpretation
            dx(2,1) = - omega*r +  (gamma + sigma*(dx1 > 0 && dx3 > 0))*w + (delta + tau*(dx1 > 0 && dx3 > 0))*y;
        case 5
            %% Situation 5: William's interpretation
            ws      = 0.01;     % transition width for sigma (units: same as sigma)
            wt      = 0.01;     % transition width for tau (units: same as tau)
            dx(2,1) = - omega*r +  (gamma + sigma/(1 + exp(-dx3/ws)))*w + (delta + tau/(1 + exp(-dx1/wt)))*y;    
        otherwise
            warning('!!!')
    end

    dx(3,1) = - theta*w + lambda*y*w - mu*w*w;
end

[Christopher]

Ok, this is incredibly helpful! I am just trying to figure out how to show that the damped oscillations can cause the wage share to eventually fall to zero. Effectively, what I am trying to do is show that the wage share falls more during recessions than it rises during booms in the cycle. It looks like from the 2d phase plane plot that the wage share cycles around some positive value. I am trying to demonstrate cycles that tend toward zero in output, wage share space.

Answer 2

[Torsten]

You mean

% linear ODEs
dx(1,1) = alpha * y - beta * r * y;    
dx(3,1) = - theta * w + lambda * y * w - mu * w * w;
if dx(1,1) > 0 && dx(3,1) > 0
    dx(2,1) = - omega * r + gamma * w + sigma * w + delta * y + tau * y;
else
    dx(2,1) = - omega * r + gamma * w  + delta * y;
end

?

[Christopher]

Thanks for this! That might work. I just want the parameters sigma and tau to react (geater than zero) to wdot and ydot when they are greater than zero. Does that make sense?

[Torsten]

I don't know if it makes physical sense.

I just technically did what you wanted: I set sigma and tau to 0 if at least one of ydot or wdot is <= 0.

[Christopher]

@Torsten Yes, not asking in a physical sense...

[Christopher]

@Torsten I am getting this error when I enter the code.

Error: File: IFStatementAdded.m Line: 59 Column: 1

All functions in a script must be closed with an 'end'.

[Torsten]

Works for me. I also included the changes to be made for @William Rose interpretation of your question.

close all
clear all
clc
tspan   = [0 100];
y0      =  0.7;   % initial value of state variable x1
r0      =  0.7;   % initial value of state variable x2
w0      =  0.8;   % initial value of state variable x3
x0      = [y0; r0; w0];
[t, x]  = ode45(@odefcn, tspan, x0);
%% Plot x2 (y-axis) vs. x1 (x-axis)
subplot(3,1,1);
plot(t,x(:,1));
xlabel('Time')
ylabel('Output');
subplot(3,1,2);
plot(t,x(:,2));
xlabel('Time')
ylabel('Policy Rate');
subplot(3,1,3);
plot(t,x(:,3));
xlabel('Time')
ylabel('Wage Share');

y0 = [0.07; 0.07; 0.8];
%% System of three linear differential equations
function dx = odefcn(t, x)
    % definitions
    y       = x(1);
    r       = x(2);
    w       = x(3);
    
    % parameters
    alpha   = 1.0;
    beta    = 1.0;
    gamma   = 0.3;
    delta   = 0.3;
    mu      = 0.4;
    lambda  = 1.0;
    theta   = 0.4;
    omega   = 0.4;
    sigma   = 0.1;
    tau     = 0.1;
    
    % linear ODEs
    dx(1,1) = alpha * y - beta * r * y;    
    dx(3,1) = - theta * w + lambda * y * w - mu * w * w;
        
    % My interpretation of your question
    if dx(1,1) > 0 && dx(3,1) > 0
      dx(2,1) = - omega * r + gamma * w + sigma * w + delta * y + tau * y;
    else
      dx(2,1) = - omega * r + gamma * w  + delta * y;
    end
    
    % William Rose's interpretation of your question
    %if dx(1,1) <= 0
    %  tau = 0;
    %end
    %if dx(3,1) <= 0
    %  sigma = 0;
    %end
    %dx(2,1) = - omega * r + gamma * w + sigma * w + delta * y + tau * y;
end

[Sam Chak]

Frankly speaking, is the system still considered as a linear system given that dx(3) has a quadratic term?

Additionally, due to the relatively small magnitudes of sigma and tau, the change from (gamma + sigma)*w to gamma*w, and the change from (delta + tau)*y to delta*y seem to have insignificant effects on the responses.

[William Rose]

Good points, @Sam Chak.

Answer 3

[Walter Roberson]

Do not use an if statement inside an ODE...

Not unless

• the if statement will always follow the same branch for any one invocation of the ode function; OR
• the second derivative of the branches is continuous

sigma > 0 , tau > 0 if dx(3,1) > 0, dx(1,1) > 0 and sigma = 0 , tau = 0

That does not have a continuous second derivative; it probably does not have a continuous first derivative either.

In your particular case, it might make sense to impose 'NonNegative' https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html -- but that depends on whether your conditions are "or" or "and"

What you should be doing is using an ODE event function to determine when the conditions are met, with the "isterminal" set to true, so that you leave the ode*() call when the conditions are true. Then modify boundary conditions as necessary and call the ode*() function again to pick up at the time you left off.

Answer 4

[William Rose]

@Christopher,

[edit: fix spelling error; no changes to code]

[Edit: Change "./" to "/" inside function odefcn(), since "./" is not needed. Should work either way.]

You say "I am trying to create an IF statement sigma > 0 , tau > 0 if dx(3,1) > 0, dx(1,1) > 0 and sigma = 0 , tau = 0 otherwise"

In the current version of your code, sigma=0.1 and tau=0.1. I assume you want sigma=0.1 when dx(3)>0, otherwise sigma=0. And I assume you want tau=0.1 when dx(1)>0, tau=0 otherwise. If that is not true, then please be more specific in your question.

I have dealt with this same issue when modeling the cardiovascular system: there are one-way valves in the system which allow forward flow, but not backflow. I learned that ode45() and other solvers do not like "if" statements.

Make sigma be a smooth (i.e. sigmoidal) function of dx(3).

where ws is the transition width for sigma*. Do the same for tau.

Here's a plot of sigma(dx(3)) when ws=0.01. I will also compute and plot sigma shifted to the right by 2*ws, in case you want the transition to happen on the positive side of zero:

ws=0.01; dx3=-1:.005:1; 
sigma=0.1./(1+exp(-dx3/ws));
sigma2=0.1./(1+exp(-(dx3-2*ws)/ws));
plot(dx3,sigma,'-b.',dx3,sigma2,'--b'); grid on; xlabel('dx_3'); ylabel('\sigma')
legend('sigma','sigma2')

Inside odefcn():

• Comment out the lines that assign constant values to sigma and tau.
• Compute dx(1) and dx(3) (they don't depent on sigma or tau), but do not compute dx(2).
• Compute sigma and tau using sigmoid functions of dx(3) and dx(1) respectively.
• Compute dx(2) (depends on sigma and tau).

It will be easier for others to assist you if you format your code as code (highlight code section, then click the "code" icon at top of window).

%% System of three linear differential equations
function dx = odefcn(t, x)
% definitions
y       = x(1);     % output
r       = x(2);     % policy rate
w       = x(3);     % wage share

% parameters
alpha   = 1.0;
beta    = 1.0;
gamma   = 0.3;
delta   = 0.3;
mu      = 0.4;
lambda  = 1.0;
theta   = 0.4;
omega   = 0.4;
% sigma   = 0.1;
% tau     = 0.1;
ws = 0.01;  % transition width for sigma (units: same as sigma)
wt = 0.01;  % transition width for tau (units: same as tau)
% linear ODEs
dx=zeros(3,1); % allocate column vector
dx(1) = alpha * y - beta * r * y;
dx(3) = - theta * w + lambda * y * w - mu * w * w;
sigma = 0.1/(1+exp(-dx(3)/ws));
tau   = 0.1/(1+exp(-dx(1)/wt));
dx(2) = - omega * r + gamma * w + sigma * w + delta * y + tau * y;
end

* Transition width, ws: In width 2*ws, sigma goes from 1/(1+e)=27% to e/(1+e)=73% of its final value. If you cannot tolerate any significant positive values of sigma when dx(3) is <=0, then you can make ws smaller or you can slide the function to the right a bit - see plot above.

[William Rose]

@Christopher,

You indicate in your comment to @Torsten that you want sigma and tau to depend on the "and" of dx(1) and dx(3). If you are still having trouble with if statements, then you can do something like

dx1=-1:.005:1; dx3=dx1; % this line would not be in your code - just to demonstrate
w=0.01;
sigma=0.1./((1+exp(-dx1'/w))*(1+exp(-dx3/w)));
surf(dx1,dx3,sigma,'EdgeColor','none')
xlabel('dx1'); ylabel('dx3'); zlabel('\sigma')

●

The plot shows that sigma=1 when dx1 and dx3 are both >0, and sigma=0 otherwise. Do likewise for tau. You can plug that idea into your odefcn().