How to apply If else condition properly for ode45 function file?

Hi, I am trying to solve system of first order ode via ode45. Please find the attached file for it. My system of equations are

x(2,1) %to show what is the value of my variable.
if x(2,1)+omega>0
  dx(1,1)=x(2,1); %These are the first set of equations of an unstable oscillator (negative   %damping).
  dx(2,1)=0.01*x(2,1)-x(1,1);
else
  keyboard 
  dx(1,1)=0;
  dx(2,1)=0;
end

Now when i m using this file in ODE45 for omega=2, then as the system of equations is changed from 1 to 2, there should be no change in the solution of x(1,1)and x(2,1), but during the running they are changing. Please help me out for this.

Regards Sunit

Answer 1

[Mischa Kim]

Sunit, use something like

 function myODE(omega)
 IC    = [1 -1];
 options = odeset('AbsTol',1e-10,'RelTol',1e-10);
 [t,X]   = ode45(@test,[0 2],IC,options,omega);
 figure
 plot(t,X(:,2),'r',t,X(:,1),'b')
 xlabel('t')
 ylabel('x(1), x(2)')
 grid
 end

 function dx = test(t,x,omega)
 if (x(2) + omega)>0
    dx = [x(2,1); ...
          0.01*x(2,1)-x(1,1)];
 else
    dx = [0; 0];
 end
 end

Is this the behavior you are looking for?

[Sunit Gupta]

Hi Mischa, Thanx for ur answer, but this is not what i m looking for. As the system of equations changed from 1 to 2, ideally the values of x(1) and x(2) should not change. But they are changing. To test it, just put

function dx = test(t,x,omega)
x(2)
 if (x(2) + omega)>0
    dx = [x(2,1); ...
          0.01*x(2,1)-x(1,1)];
 else
    keyboard
    dx = [0; 0];
 end
 end

Now as the system will change from 1 to 2 you can see that the values of variables will change.

[Mischa Kim]

I am not sure I understand. Let's say you start with initial conditions such that you end up with system 1 (the non-trivial ODE). When x(2) + omega goes negative, the code switches to system 2. From that time on the state vector remains constant. With the above IC, use e.g.

 myODE(1.2)

and remove the semi-colons from the two dx = ... commands. If you change the solver settings and replace the ode call by

 options = odeset('AbsTol',1e-10,'RelTol',1e-10);
 [t,X]   = ode45(@test,[0 10],IC,options,omega);

you'll see that the state vector remains constant when x(2) reaches -1.2.

On the other hand if you execute

 myODE(1)

the ode solver starts with system 2 (as it should) and remains in that state (as it should).

Answer 2

[Parham]

I tested this but the thing you're saying did not happen:

function []=pap()

[AA,BB] = ode45(@testttt,[0,5],[2,.1])
end
function dx= testttt(t,x)
x(2,1)
omega = 2;
if x(2,1)+omega>0
    dx(1,1)=x(2,1);
    dx(2,1)=0.01*x(2,1)-x(1,1);
else
    dx(1,1)=0;
    dx(2,1)=0;
end

end

[Sunit Gupta]

Hi, That did not happen as the time span is too short. Make it run for [0 300], then probably you can see.

Regards Sunit

[Mischa Kim]

See plot screen shot above for myODE(1.2) and IC = [1 -1].

[Parham]

Ok. The problem is when MATLAB integrates, it finds the jacobian and use it while integrating. That is the reason why the time increments of the beginning of the integration is very small (0.01) and increase to larger values ( 10 20 ...). So, when the integration time increases, your value (x(1,2)) might go below the omega value that you defined and not exactly stop at that value. This kind of switching I guess can be handled using the event handling of the ode solver. Check out the event handling of ode solvers.

[Sunit Gupta]

Thanx Parham, I think the same. It is happening because of numerical integration only. But can i handle this by fixing the time step, say 0.001?

Regards Sunit