Using event function to switch between two sets of ODEs

Hi,

I wrote a code to solve two sets of ODEs, one valid for r < a and the other valid for r≥a. I used the event function to switch between the two. However, it is not triggering. Any help would be great. Here is the code:

clear all;

R = 0.5;
St =0.01;
A = R/St;             
U = 1;
a = 1;

tspan = [0 5.5];
tstart = tspan(1);
tend = tspan(end);

y0 = [0.1 0.1 0.1 0.1 0.1 0.1];

% Accumulators for ODE solve output
tout = tstart;  % global solution time
yout = y0;      % global solution position
teout = [];     % time at which events occur
yeout = [];     % position at which events occur
ieout = [];     % flags which event trigger the switch

options = odeset('Events',@(t,y)RV_event1(t,y,a));
% Main loop
while tout(end) < tend
     if  yout(end,1)^2+yout(end,3)^2+yout(end,5)^2 < a 
         disp('hello')
         [t,y,te,ye,ie] = ode15s(@(t,y)HS_car_in(t,y,A,R,U,a),tspan ,y0,options);  
     else
         disp('goodbye')
         [t,y,te,ye,ie] = ode15s(@(t,y)HS_car_ou(t,y,A,R,U,a),tspan ,y0,options);
     end
     nt = length(t);
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,:)];
    teout = [teout; te]; % Events at tstart are never reported.
    yeout = [yeout; ye];
    ieout = [ieout; ie];
    tspan = [t(end), tend]; % reset time span,'where you just ended to tend'
    y0 = [y(end,1);y(end,2);y(end,3);y(end,4);y(end,5);y(end,6)];           % reset IC as where you just ended
end

[phi,theta,r] = cart2sph(y(:,1),y(:,3),y(:,5));
polarplot(theta,r,'k')

function du = HS_car_in(t,u,A,R,U,a)
du = zeros(6,1);
du(1) = u(2);
du(2) =  -A*u(2) - 3*((9*R*(u(1)^2 + u(3)^2 - a^2/2)*U)/2 + A*u(5)*a^2)*u(1)*U/(2*a^4);
du(3) = u(4);
du(4) =  -A*u(4) - 3*((9*R*(u(1)^2 + u(3)^2 - a^2/2)*U)/2 + A*u(5)*a^2)*U*u(3)/(2*a^4);
du(5) = u(6);
du(6) = -A*u(6) + 3*U*(-A*a^4/2 + ((u(1)^2 + u(3)^2 + u(5)^2/2)*A - (9*R*U*u(5))/4)*a^2 + (9*R*U*u(5)^3)/4)/a^4;
end

function du = HS_car_ou(t,u,A,R,U,a)
du = zeros(6,1);
du(1) = u(2);
du(2) = -A*u(2) -(3*u(1)*((A*u(5)^3 - 6*U*u(5)^2*R + A*(u(1)^2 + u(3)^2)*u(5) + (3*U*R*(u(1)^2 + u(3)^2))/2)*(u(1)^2 + u(3)^2 + u(5)^2)^(3/2) + (3*R*U*a^3*(u(1)^2 + u(3)^2 + 5*u(5)^2))/4)*U*a^3)/(2*(u(1)^2 + u(3)^2 + u(5)^2)^5); 
du(3) = u(4);
du(4) = -A*u(4) - 3*((A*u(5)^3 - 6*U*u(5)^2*R + A*(u(1)^2 + u(3)^2)*u(5) + (3*U*R*(u(1)^2 + u(3)^2))/2)*(u(1)^2 + u(3)^2 + u(5)^2)^(3/2) + (3*R*U*a^3*(u(1)^2 + u(3)^2 + 5*u(5)^2))/4)*U*u(3)*a^3/(2*(u(1)^2 + u(3)^2 + u(5)^2)^5);
du(5) = u(6);
du(6) = -A*u(6) + U*((u(1)^2 + u(3)^2 + u(5)^2)^(3/2)*((u(1)^2 + u(3)^2 + u(5)^2)*(u(1)^2 + u(3)^2 - 2*u(5)^2)*A - (27*R*U*u(5)*(u(1)^2 + u(3)^2 - (2*u(5)^2)/3))/2)*a^3/2 + (u(1)^2 + u(3)^2 + u(5)^2)^5*A - (9*R*U*a^6*u(5)^3)/2)/(u(1)^2 + u(3)^2 + u(5)^2)^5;
end

function [value,isterminal,direction] = RV_event1(t,y,a)
value = sqrt(y(1)^2+y(3)^2+y(5)^2) - a ;     % crossing condition, evaluates to zero at event condition
isterminal = 0;                                      
direction = 0;
end

[Walter Roberson]

isterminal = 0;                                      

The event needs to be terminal.

[Takumi]

The conditions appear to be different, is this intentional?

yout(end,1)^2+yout(end,3)^2+yout(end,5)^2 < a

sqrt(y(1)^2+y(3)^2+y(5)^2) - a

[Walter Roberson]

Ah, they differ in sqrt() or not.

[Nassima]

Walter,

When I use isterminal = 1, it only solves the second set of ODEs. What I want is to solve the first set of ODEs while r<a. Then as soon as the solution reaches r=a, I want the code to switch to the second set of ODEs i.e. r≥a.

As for the above conditions, a=1 and so the sqrt does not have an effect. Is it what you mean?

Thanks.

[Walter Roberson]

When I test with your code, it runs the first one ('hello') until an appropriate event where there is a zero crossing on the constraint (that is, r becomes equal to a).

Then it runs the second ode ('goodbye'). It is starting from the point r == a (to within roundoff) because that was the termination condition for the first series. When you start an ode with an event, it records the initial value of the constraint, and it does not trigger again until the constraint value changes -- not until there is another crossing. But again the termination is r == a. But your condition to run the first ode is

yout(end,1)^2+yout(end,3)^2+yout(end,5)^2 < a

with a strict inequality, and so when that expression == a instead of < a then the first ode is not run, so it fires up the second ode. And that is going to continue to happen the same way, that the second ode continues to get invoked, because you did not program alternating them, you programmed a strict inequality and a termination condition that searches internally hoping to find strict equality, only giving up and getting a non-zero value from the subtraction if the step size falls below the tolerance.

ODE events do not terminate the ode as soon as they first find a location that causes a sign change in the constraint: ODE events have Regula-Falsi search code to try to find the location at the zero crossing, where ideally sqrt(y(1)^2+y(3)^2+y(5)^2) == a exactly. They then terminate at that time and boundary conditions. If you want to know which direction you were going, you cannot figure that out just from the last output boundary values.

What you could potentially do is look at the second-last output boundary values: if yout(end-1,1)^2+yout(end-1,3)^2+yout(end-1,5)^2 < a then you must have gone up from a lower value.

Or you could just toggle between the cases, unless you are expecting the process to just "touch" a and then reverse direction.

Answer 1

[Takumi]

How about conditional switching in a differential equation?

clear all;
R = 0.5;
St =0.01;
A = R/St;             
U = 1;
a = 1;
tspan = [0 5.5];
y0 = [0.1 0.1 0.1 0.1 0.1 0.1];

 [t,y] = ode15s(@(t,y)HS_car(t,y,A,R,U,a),tspan ,y0);  
[phi,theta,r] = cart2sph(y(:,1),y(:,3),y(:,5));
polarplot(theta,r,'k')

function du = HS_car(t,u,A,R,U,a)
if sqrt(u(1)^2+u(3)^2+u(5)^2)<a     % or y(1)^2+y(3)^2+y(5)^2 <a ???
    du(1) = u(2);
    du(2) =  -A*u(2) - 3*((9*R*(u(1)^2 + u(3)^2 - a^2/2)*U)/2 + A*u(5)*a^2)*u(1)*U/(2*a^4);
    du(3) = u(4);
    du(4) =  -A*u(4) - 3*((9*R*(u(1)^2 + u(3)^2 - a^2/2)*U)/2 + A*u(5)*a^2)*U*u(3)/(2*a^4);
    du(5) = u(6);
    du(6) = -A*u(6) + 3*U*(-A*a^4/2 + ((u(1)^2 + u(3)^2 + u(5)^2/2)*A - (9*R*U*u(5))/4)*a^2 + (9*R*U*u(5)^3)/4)/a^4;
else
    du(1) = u(2);
    du(2) = -A*u(2) -(3*u(1)*((A*u(5)^3 - 6*U*u(5)^2*R + A*(u(1)^2 + u(3)^2)*u(5) + (3*U*R*(u(1)^2 + u(3)^2))/2)*(u(1)^2 + u(3)^2 + u(5)^2)^(3/2) + (3*R*U*a^3*(u(1)^2 + u(3)^2 + 5*u(5)^2))/4)*U*a^3)/(2*(u(1)^2 + u(3)^2 + u(5)^2)^5); 
    du(3) = u(4);
    du(4) = -A*u(4) - 3*((A*u(5)^3 - 6*U*u(5)^2*R + A*(u(1)^2 + u(3)^2)*u(5) + (3*U*R*(u(1)^2 + u(3)^2))/2)*(u(1)^2 + u(3)^2 + u(5)^2)^(3/2) + (3*R*U*a^3*(u(1)^2 + u(3)^2 + 5*u(5)^2))/4)*U*u(3)*a^3/(2*(u(1)^2 + u(3)^2 + u(5)^2)^5);
    du(5) = u(6);
    du(6) = -A*u(6) + U*((u(1)^2 + u(3)^2 + u(5)^2)^(3/2)*((u(1)^2 + u(3)^2 + u(5)^2)*(u(1)^2 + u(3)^2 - 2*u(5)^2)*A - (27*R*U*u(5)*(u(1)^2 + u(3)^2 - (2*u(5)^2)/3))/2)*a^3/2 + (u(1)^2 + u(3)^2 + u(5)^2)^5*A - (9*R*U*a^6*u(5)^3)/2)/(u(1)^2 + u(3)^2 + u(5)^2)^5;
end
du = du';
end

[Walter Roberson]

No, you should not do this. Event functions are the way to go. The ode* routines require that the functions are continuously differentiable twice, and if that is not true then you will get incorrect answers (or it will complain about being unable to integrate.)

[Takumi]

Thank you for commenting and I'm sorry for answering incorrectly.

[Nassima]

Thank you Takumi. I have already tried it and Walter is right, event function is the only way.