How do I stop ode solver when output is complex number?

15 次查看（过去 30 天）
Vikash Pandey
Vikash Pandey 2018-11-15
评论： Umar Abdurrehman 2021-2-16
I am currently working on a coupled ode problem, and using ode45 solver to solve it.
My code is something like this:
bubble.m
function rdot = f(t, r)
%Some mathematical expressions
....
...
...
rdot(1) = r(2); % first derivative of r1
rdot(2) = (X11 - X21 - (X31*(X41 - X51)))/X81; % second derivative of r1
rdot(3) = r(4); % first derivative of r2
rdot(4) = (X12 - X22 - (X32*(X42 - X52)))/X82; % second derivative of r2
rdot = rdot';
And,
bubble_plotter.m
clc;
clear all;
%close all;
time_range = [0 1d-4];
r1_eq = 10d-6;
r2_eq = 1d-6;
initial_conditions = [r1_eq 0.d0 r2_eq 0.d0];
[t, r] = ode45('bubble', time_range, initial_conditions);
r1_us=1000000*r(:,1);
r1dot=r(:,2);
r2_us=1000000*r(:,3);
r2dot=r(:,4);
%% Recovery of double derivatives
for i=1:numel(t)
t_actual = t(i);
r_actual = r(i,:);
rdot(i,:) = bubble_mettin(t_actual,r_actual);
end
r1ddot = rdot(:,2);
r2ddot = rdot(:,4);
figure(101);
plot(normalized_time, r1_us, 'b-', normalized_time, r2_us, 'k-.');
For certain input parameter values, ode solver successfully runs only for first few t-values and outputs corresponding real positive values for r, which in the code is r(1) and r(3). But after that the solver starts outputing negative or complex values (that is, imaginary numbers) for r(1) and r(3). I wish the program stops right before it starts yielding complex or negative values. Also the so far calculated values of r(1), r(2), r(3), r(4), rdot(2) and rdot(4) must be receoverable.
Another condition to stop the solver could be that r(1) + r(3) should not excced a given value.
  4 个评论
显示 2更早的评论

请先登录，再进行评论。

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

回答（1 个）

madhan ravi
madhan ravi 2018-11-15
编辑：madhan ravi 2018-11-15
  5 个评论
显示 3更早的评论
madhan ravi
madhan ravi 2018-11-15
Opt = odeset('Events', @myEvent); %add this before ode45
[t, r] = ode45('bubble', time_range, initial_conditions);
save the below function separately:
function [value, isterminal, direction] = myEvent(t, r)
value = isreal(r(1))==0 | r(1)<0 | isreal(r(3))==0 | r(3)<0;
isterminal = 1; % Stop the integration
direction = 0;
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by