I was searching for ways to terminate MATLAB ode when it meets a particular condition. I found the answer in this topic MatLab ODE start/stop conditions, where the use of 'events' was discussed. However this applies to ode45 and when I attempted to use 'events' with ode15i, it simply didn't worked and MATLAB is showing error.
I am trying to learn this with simple example and solved a simple system of differential equations as given below.
dx/dt = 5x + 3y; dy/dt = x + 7y; I solved them using ode45 and tried to do the same with ode15i, but it doesn't work. Given below are my codes.
With ode45
function start_stop_test_ode45
y0 = [5;1];
tv = linspace(0,2,100);
options = odeset('Events',@events);
f = @(t,y) [5*y(1) + 3*y(2);y(1) + 7*y(2)];
[t,Y] = ode45(f,tv,y0,options);
xNI = Y(:,1);
yNI = Y(:,2);
xCF = 3*exp(4*t) + 2*exp(8*t);
yCF = -1*exp(4*t) + 2*exp(8*t);
% Here we plot all the graphs
figure(1)
plot(t,xNI,'--k',t,xCF,'r','Linewidth',1.75)
xlabel('t (s)')
ylabel('x')
legend('Numerical Solution','Closed Form Solution')
figure(2)
plot(t,yNI,'--k',t,yCF,'r','Linewidth',1.75)
xlabel('t (s)')
ylabel('y')
legend('Numerical SOlution','Closed Form Solution')
% Here we solve plot the variation of x with y
figure(3)
plot(xNI,yNI,'k','Linewidth',2);
end
function [value,isterminal,direction] = events(t,y)
value = [y(1) - 7782;y(2) - 8863]; % Detect y = 7356
isterminal = [1;1];
direction = [0;0];
end
With ode15i
function start_stop_test_ode15i
clc;clear all
t0 = 0;
y0 = [5;1];
Fxdy0 = [1;1];
Fxdyp0 = [0;0];
yp0 = [28;12];
tRange = [0 2];
options = odeset('Events',@events);
[y0,yp0] = decic(@ode15ifun,t0,y0,Fxdy0,yp0,Fxdyp0);
sol = ode15i(@ode15ifun,tRange,y0,yp0,options);
tv = linspace(0,2,100);
sv = deval(sol,tv);
sv = sv';
t = tv;
xNI = sv(:,1);
yNI = sv(:,2);
xCF = 3*exp(4*t) + 2*exp(8*t);
yCF = -1*exp(4*t) + 2*exp(8*t);
% Here we plot all the graphs
figure(4)
plot(t,xNI,'--k',t,xCF,'r','Linewidth',1.75)
xlabel('t (s)')
ylabel('x')
legend('Numerical Solution','Closed Form Solution')
figure(5)
plot(t,yNI,'--k',t,yCF,'r','Linewidth',1.75)
xlabel('t (s)')
ylabel('y')
legend('Numerical SOlution','Closed Form Solution')
% Here we solve plot the variation of x with y
figure(6)
plot(xNI,yNI,'k','Linewidth',2);
end
function [value,isterminal,direction] = events(t,y)
value = [y(1) - 7782;y(2) - 8863]; % Detect y = 7356
isterminal = [1;1];
direction = [0;0];
end
Where the ode15ifun is
function res = ode15ifun(t,y,yp)
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
res1 = yp(1) - 5*y(1)- 3*y(2);
res2 = yp(2) - y(1) - 7*y(2);
res = [res1;res2];
end
ode45 is working fine but while using ode15i I am getting error message. Can anyone help as to how to do the same with ode15i?
Thank you very much
(In Response to TroyHaskin) I am adding an image file of the error message for your reference.
