I am trying to incorporate multiple IF statements in my ODE to generate a single output.

There are four cases that I am trying to incorporate as multiple IF statements. Here are the four cases below. The first case is the one I used in my code to generate the output.

dx(1,1) > 0 && dx(3,1) > 0

dx(1,1) > 0 || dx(3,1) > 0

dx(1,1) > 0 && dx(3,1) <= 0

dx(1,1) <= 0 && dx(3,1) > 0

Below is the code that I used to generate the first case:

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 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');

y0 = [0.7; 0.7; 0.8];

%% System of three 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.6;
    delta   = 0.6;
    mu      = 0.1;
    lambda  = 0.6;
    theta   = 0.4;
    omega   = 0.4;
    sigma   = 0.1;
    tau     = 0.1;
    
    % ODEs
    dx(1,1) = alpha * y - beta * r * y;    
    dx(3,1) = - theta * w + lambda * y * w - mu * w * w;
    
    % Asymmetrical Reaction Function
    if dx(1,1) > 0 && dx(3,1) > 0
        dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r;
    else
        dx(2,1) = - omega * r + gamma * w * r + delta * y * r;
    end
end

[Torsten]

And what is your question ?

[Christopher]

How do I code multiple IF statements to represent the three remaining cases?

[Christopher]

@Sam Chak Thanks for adding the output!

[Sam Chak]

Hi @Christopher

Please try updating your code using this structure and see if it functions as expected. Ideally, you should ensure continuities between the cases, when running the ode45 solver.

if     case1
    run this
elseif case2
    run this
elseif case3
    run this
elseif case4
    run this
end

[Torsten]

dx(1,1) = 5;
dx(3,1) = -2;
if dx(1,1) > 0
    if dx(3,1) > 0
        c = 'case 1'
    else
        c = 'case 3'
    end
else
    if dx(3,1) > 0
        c = 'case 4'
    else
        c = 'not included in your cases: dx(1,1) <= 0 and dx(3,1) <= 0'
    end
end

c = 'case 3'

[Torsten]

@Sam Chak

Case 1 and Case 2 are not mutually exclusive - this will cause different results depending on the order in which the if-statements are processed.

[Sam Chak]

@Torsten, you are correct. I did not assess the logical conditions previously.

[Christopher]

@Sam Chak @Torsten This is helpful. Thank you both! I was able to run something smilar to what @Sam Chak had. I understand @Torsten's point, but I am not sure how exactly to incorporate the code into my exiting ODE's. Here is what I did which I understand is not accurate:

% ODEs
dx(1,1) = alpha * y - beta * r * y;    
dx(3,1) = - theta * w + lambda * y * w - mu * w * w;

% Asymmetrical Reaction Function
if     dx(1,1) >  0 && dx(3,1) >  0
    dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r;
elseif dx(1,1) >  0 || dx(3,1) >  0
    dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r;
elseif dx(1,1) >  0 && dx(3,1) <= 0
    dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r;
elseif dx(1,1) <= 0 && dx(3,1) >  0
    dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r;
else
    dx(2,1) = - omega * r + gamma * w * r + delta * y * r;
end

[Walter Roberson]

All of the branches are the same except for the last one. You could compact the if tree to

if dx(1,1) <= 0 && dx(3,1) <= 0
    dx(2,1) = - omega * r + gamma * w * r + delta * y * r;
else
   dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r; 
end

Answer 1

[Sam Chak]

Hi @Christopher

Based on your definitions of the four cases in this comment and @Walter Roberson's input, we can run the simulation to observe the effect of the discontinuous term, , on the 2nd differential equation for the asymmetrical reaction.

tspan   = linspace(0, 100, 100001);
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 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');

dx  = zeros(numel(t), 3);
term  = zeros(numel(t), 1);
for i = 1:numel(t)
    [dx(i,:), term(i,:)]  = odefcn(t', x(i,:)');
end

%% Observe the effect of discontinuous term
figure
subplot(211)
plot(t, term),    grid on, title({'Effect of discontinuous term, $\delta y r + \tau y r$'}, 'interpreter', 'latex', 'fontsize', 12)
ylim([-1, 5])
subplot(212)
plot(t, dx(:,2)), grid on, title('dx_2')
ylim([-2, 4])

●

%% System of three differential equations
function [dx, term] = odefcn(t, x)
    % definitions
    y       = x(1);
    r       = x(2);
    w       = x(3);
    
    % parameters
    alpha   = 1.0;
    beta    = 1.0;
    gamma   = 0.6;
    delta   = 0.6;
    mu      = 0.1;
    lambda  = 0.6;
    theta   = 0.4;
    omega   = 0.4;
    sigma   = 0.1;
    tau     = 0.1;
    
    % 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
        term    = 0;
    else
        term    = delta * y * r + tau * y * r;
    end
    
    % for comparison without if–else
    % dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r;
    
    % Asymmetrical Reaction Function
    if dx(1,1) <= 0 && dx(3,1) <= 0
        dx(2,1) = - omega * r + gamma * w * r + delta * y * r;
    else
      % dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + delta * y * r + tau * y * r; 
        dx(2,1) = - omega * r + gamma * w * r + sigma * w * r + term; 
    end
end

[Christopher]

Ok, this makes perfect senese! Is it possible to plot all three equations together? I just want to see if the cyclical behavior is still intact with the discontinuity?

[Christopher]

@Sam Chak Nevermind! It looks like you have this plotted already!

[Christopher]

@Sam Chak @Walter Roberson @Torsten This looks great!

Answer 2

[Walter Roberson]

The mathematics of ode45 and similar routines is such that the second derivative of the input expressions must be continuous over any one call to the ode*() routine.

If the first derivative is not continuous, then about 1/3 of the time, ode45 detects that and issues an error message.

If the second derivative is not continuous, then typically ode45 does not notice and simply gives the wrong answer.

It is rather uncommon for ode routines that contain if statements to be continuous in the second derivative. The sample code you provided is not continuous in the second derivative.

The right way to handle this situation is to use ode event functions to determine the transition boundaries, marking each transition as "terminal", and to loop over ode45() calls, recording the partial outputs produced, and calling ode45() again to pick up from where the previous call left off.