2

I have a system of ODE equations which I want to solve, but there is a tricky part that when the system reaches steady-state, I would like to change the value of one (or more) parameters. For example, consider the following:

function dydt = diff(t,x,params)
   F = params(1);
   G = params(2);
   dydt = zeros(2,1);
   dydt(1) = F*x(1) - G*x(1)*x(2);
   dydt(2) = (F-G)*x(2);
end

I would like my code to work such that when the system has reached steady-state, the value of F is changed to 10 and the value of G is changed to 2, for example. I was thinking of detecting the values of dydt(1) and dydt(2) by using, for example, 

if norm(dydt)<1
   F = 10;
   G = 2;
end

How do I do that for ODE expression in Matlab? If I put this if condition before the ODE expression, I the value of dydt will always be zero. But if I put this If condition after the ODE expression, the If conditions will have no use to correct the ODE expression.

Thank you!

user6101792
  • 53
  • 7
  • 1
    The proper way to do this is to use [event location](http://www.mathworks.com/help/matlab/math/ode-event-location.html) to stop the integration based on your steady-state condition ([example here](http://stackoverflow.com/a/15992197/2278029)). Then change your parameters and start a new integration. Avoid putting an `if` statement in your ODE function to change parameters. – horchler Apr 17 '16 at 15:39
  • I considered event location, but does this mean that I have to stop the integration and then start over a new integration procedure? If so, do I have to update the initial values to the last points of the first integration step before steady-state was found? – user6101792 Apr 18 '16 at 02:49

1 Answers1

0

A parameter of an ODE is assumed to be fixed and does not depend on the state of the system. What you're trying to do is simulate a piecewise continuous ODE. This sort of problem is usually solved with event location (assuming a solution exists). You need to stop the simulation at the key point, change your parameters, and restart a new simulation with initial conditions the same as those at the end of the previous.

Here's a little example with your ODE function. I don't know your initial conditions, initial parameter values, or other settings, so this is just to demonstrate this scheme:

function eventsdemo
params = [-1.5 1];
tspan = [0 10];
x0 = [1;1];
opts = odeset('Events',@events);
[t1,x1] = ode45(@(t,x)f(t,x,params),tspan,x0,opts); % Simulate with events

% Change parameters, set initial conditions based on end of previous run
params = [1.5 1];
x0 = x1(end,:);
tspan = [t1(end) 10];
[t2,x2] = ode45(@(t,x)f(t,x,params),tspan,x0); % Simulate again

% Concatenate results, removing duplicate points
t = [t1;t2(2:end)];
x = [x1;x2(2:end,:)];

figure;
plot(t,x);
hold on;
plot(t2(1),x2(1,:),'k*'); % Plot event location


function dxdt = f(t,x,params) %#ok<INUSL>
F = params(1);
G = params(2);
dxdt = [F*x(1) - G*x(1)*x(2);
        (F-G)*x(2)];


function [value,isterminal,direction] = events(t,x) %#ok<INUSL>
value = norm(x)-1e-3; % Don't try to detect exact zero for asymptotics
isterminal = true;
direction = -1;

In this case the solution asymptotically approaches a stable fixed point at (0,0). It's important to not try to detect (0,0) exactly as the solution may never reach that point. You can instead use a small tolerance. However, depending on your system, it's possible that choosing this tolerance could impact the behavior after you change parameters. You could also consider reformulating this problem as a boundary value problem. I don't know what you're trying to do with this system so I can't say much else (and it would probably be off-topic for this site).

horchler
  • 18,384
  • 4
  • 37
  • 73