0

Bit of a novice here... I'm trying to a solve a system of 6 odes but with a condition that if (x(1)<=0.1) then x(6)=x(6)/2.0. How do I put that into the code? Thanks!

Fred
  • 9
  • Please state the condition in a mathematical form. x(6)=x(6)/2.0 makes sense as an assignment, but not as a condition. –  Jul 04 '15 at 02:29
  • And what if `x(1)>0.1`? What have you tried so far? From what I can tell, it sounds like you should look into [events](http://mathworks.com/help/matlab/ref/odeset.html#f92-1017470). You want to integrate to the point where the change occurs, stop, then restart a new integration with different parameters. Examples: [1](http://stackoverflow.com/q/22818556/2278029), [2](http://stackoverflow.com/q/13755352/2278029). – horchler Jul 04 '15 at 15:49
  • Basically the six odes are solved as normal if x(1)>0.1 but when it is <=0.1 then the 6th variable is halved. I am trying to recreate something that's already been done and can see from their plots that the 1st variable crosses the threshold several times during the simulation. Does this still sound like an event (can multiple events be handled)? So far I was trying to put an if statement inside the ode function... Thanks so much for your help!! – Fred Jul 05 '15 at 04:53
  • Looks like you will have two solutions, one for when `x(1) <= 0.1` where `x(6)` is halved, and the other when `x(1) > 0.1`, with `x(6)` as given. You solve for each region separately. – buzjwa Jul 05 '15 at 09:15

1 Answers1

0

This is slightly modified basic Example 1 from http://www.mathworks.com/help/matlab/ref/ode45.html

First make a new Matlab function and call it rigid.m. You can put any code inside, but try this:

function dy = rigid(t,y)
    dy = zeros(6,1);    % a column vector
    dy(1) = y(2) * y(3);
    dy(2) = -y(1) * y(3);   
    dy(3) = -0.51 * y(1) * y(2);

    dy(4) = y(5) * y(6);
    dy(5) = -y(4) * y(6);   
    dy(6) = -0.51 * y(4) * y(5);

    if(dy(6)<0.5)
        dy(6)=dy(6)/2;
    end
end

Now run these three lines:

options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[T,Y] = ode45(@rigid,[0 12],[0 1 1,0 1 1],options);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.',T,Y(:,4),'-',T,Y(:,5),'-.',T,Y(:,6),'.');

And it is done. Matlab solvers should handle quite well the discontinuity in the derivative, but it will depend on your problem. Anyway in this case if dy(6)<0.5 then it is halved.

Piotr Kusmierczyk
  • 56
  • 2