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

1 次查看（过去 30 天）
Christopher
Christopher 2025-5-9
评论： Christopher 2025-5-10
Ran in:
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
  9 个评论
显示 7更早的评论
Christopher
Christopher 2025-5-10
编辑：Sam Chak 2025-5-10
@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
Walter Roberson 2025-5-10
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

请先登录，再进行评论。

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

采纳的回答

Sam Chak
Sam Chak 2025-5-10
Ran in:
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
  3 个评论
显示 1更早的评论

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2025-5-9
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.

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by