Using IF condition with ODE
信息
此问题已关闭。 请重新打开它进行编辑或回答。
显示 更早的评论
I would need some help to check on the proper utilization of conditions inside ODE equations.
Below can be found the code. Graphes were managed to be plotted but they were not as I expected them to be.
----------------------------------------------------------------------------------------------------------------------------------------------------------------
%Constants
Iapp=40*60;
K(1)=0.001;
K(2)=10e-4;
K(3)=1;
K(4)=5;
K(5)=0.8;
% Initial Conditions
y0 = [3.949270942 7e-4 0.000281838];
% Time inverval
t_exp = [0 30 60 120 180 240 300];
% Solver;
options = odeset('AbsTol',1e-6);
[t,y] = ode15s(@manu,t_exp,y0,options,Iapp,K);
% Graphics:
figure(1)
%Figure 1-Ca
subplot(311)
plot(t,y(:,1),'r')
legend('Ca2+ theo')
title('Ca2+')
xlabel('treatment time / min')
ylabel('Concentration / M')
hold on
grid
%Figure 2-CO3
subplot(312)
plot(t,y(:,2),'r')
legend('CO32- theo')
title('CO32-')
xlabel('treatment time / min')
ylabel('Concentration / M')
hold on
grid
%Figure 3-OH
subplot(313)
plot(t,y(:,3),'r')
legend('OH- theo')
title('OH-')
xlabel('treatment time / min')
ylabel('Concentration / M')
hold on
grid
function dydt = manu(t,y,Iapp,K)
if y(2)>=K(2)
Ca = -K(1)*y(1);
CO3 = -K(4)*y(2);
OH = 0;
elseif y(2)<=K(2)
Ca = 0;
CO3 = K(3);
OH = Iapp/96500*K(5);
end
dydt = [Ca; CO3; OH];
end
9 个评论
darova
2020-3-30
How should check this? Where are original equations?
Faidzul Hakim
2020-3-30
darova
2020-3-30
darova
2020-3-30
Faidzul Hakim
2020-3-30
Torsten
2020-3-30
As written, your integration will never be successful. Once you reach condition 1 from condition 2, y(2) starts to decrease again and condition 2 again becomes active. So you will switch from one condition to the other over and over again.
Maybe you meant that once condition 1 is met, you work with the equations for this condition for the rest of the time ?
Faidzul Hakim
2020-3-30
Torsten
2020-3-30
In this case Walter's suggestion to use Matlab's event option is the way to go.
Or - if the equations remain that simple - you can integrate using pencil and paper.
Faidzul Hakim
2020-3-30
回答（2 个）
Walter Roberson
2020-3-30
0 个投票
You need to use Events in the options to terminate the ode integration, and then restart from that time. None of the ode* functions can deal with discontinuities in any of the derivatives (and not in two more derivatives beyond either.) The ode* functions do not always notice the discontinuities, but their mathematical models rely upon there not being any discontinuities, so even if they do not complain then the results are invalid.
4 个评论
Faidzul Hakim
2020-3-30
Walter Roberson
2020-3-30
I suggest you read through ballode
Faidzul Hakim
2020-3-30
Faidzul Hakim
2020-3-31
编辑：Faidzul Hakim
2020-3-31
Hi again Walter,
I've tried using the Events options, I seem to be stucked where I would like to redefine the initial values for all my three ODEs. Could you please show me how I can extract my value of dy(1)/dt, dy(2)/dt and dy(3)/dt when my condition of dy(2)/dt achieves 'zero'?
The illustration is given in photo below and the code can be found further down:
K(1)=0.001;
K(2)=8.5079e-4;
K(3)=1;
K(4)=5;
K(5)=0.02;
% Initial Conditions
y0 = [3.949270942 7e-4 0.000281838];
% Experimental values
t_exp = [0 30 60 120 180 240 300];
% Solver;
options = odeset('AbsTol',1e-6,'Events',@criticEvents);
[t,y,te,ye,ie] = ode15s(@manu1,t_exp,y0,options,K);
if (y(:,2)-K(2))<0
[t,y,te,ye,ie] = ode15s(@manu1,t_exp,y0,options,K);
elseif (y(:,2)-K(2))>=0
y0=[ye(:,1) ye(:,2) ye(:,3)];
[t,y,te,ye,ie] = ode15s(@manu2,t_exp,y0,options,K);
end
% Graphics:
figure(1)
%Figure 1-M
subplot(311)
plot(t,y(:,1),'r')
legend('M theo')
title('M')
xlabel('time')
ylabel('unit')
hold on
grid
%Figure 2-N
subplot(312)
plot(t,y(:,2),'r')
legend('N theo')
title('N')
xlabel('time')
ylabel('unit')
hold on
grid
%Figure 3-P
subplot(313)
plot(t,y(:,3),'r')
legend('P theo')
title('P')
xlabel('time')
ylabel('unit')
hold on
grid
function dydt = manu2(t,y,K)
M = -K(1)*y(1);
N = -K(4)*y(2);
P = 0;
dydt = [M; N; P];
end
function dydt = manu1(t,y,K)
M = 0;
N = K(3);
P = K(5);
dydt = [M; N; P];
end
function [critic,isterminal,direction] = criticEvents(t,y)
critic=(y(:,2)-K(2));
isterminal=1;
direction=+1;
end
darova
2020-4-1
Maybe you don't need event function for this case. I just add persistent variable to your ode function
function main
clear functions
% your main code
end
function dydt = manu(t,y,Iapp,K)
persistent flag
if isempty(flag)
flag = false;
end
if y(2)>=K(2) || flag
% variables
flag = true;
elseif y(2)<=K(2)
Ca = 0;
CO3 = K(3);
OH = Iapp/96500*K(5);
end
dydt = [Ca; CO3; OH];
end
result
1 个评论
Faidzul Hakim
2020-4-1
此问题已关闭。
2020-3-30
关闭：
2021-8-20
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
