MATLAB: Steady state criterion

steady stateurgent

Hi, I have a system of 4 ODE's and I solved it by using ode15s solver. But I want to stop simulations when the steady state is reached. Can someone tell me how to do it? I need to know this urgent!! Thanks and regards.

Best Answer

By "steady state" do you mean an equilibrium solution, or some non-equilibrium state that the solution settles to after an initial transient (eg a periodic solution)? If the latter, I don't know off the top of my head how you'd do that. If the former, you can use an event.

Write an event function, where, in your case the event should be something like the norm of the derivatives is zero (or, realistically, close to zero).

function [x,isterm,dir] = eventfun(t,y)
dy = copy-n-paste-from-your-ode-equation-function;
x = norm(dy) - 1e-5;
isterm = 1;
dir = 0;  %or -1, doesn't matter

This defines an event for the integration to be whenever the norm of the derivatives is less than 10^-5. You can tweak that definition according to your knowledge of the system (and the likely values of the derivatives, in particular). The value of isterm will tell ode45 to stop integrating once the event occurs. The direction shouldn't matter, because there's really no way for the norm of the derivatives to increase to 1e-5.

Once you've done this, add the event to your ode options:

opts = odeset('Events',@eventfun);

Then call ode15s with that option structure. You can use an infinite integration time (or just very big), and the integration will stop when (if!) the event occurs:

[t,z] = ode15s(@odefun,[0 Inf],y0,opts);

Related Solutions

MATLAB: Variable coefficients-help needed please!!!

I redefined your code as:

function dy  =  funl12(t,y)
dy(1,1)  =  (k1(t).*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2(t).*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;
dy(2,1)  =  (mu(t).*(y(2).^n)/(K(t)^n+y(2).^n)).*exp((-y(1)).*k)-(k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1(t)+gamma(t).*y(4)).*y(2);
dy(3,1)  =  (k3(t).*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2(t)+gamma(t).*y(4)).*y(3);
dy(4,1)  =  r(t).*y(4).*(1-(y(4)./(alpha(t).*(y(2)+y(3)-1e-12))));
      % Nested k1
      function y  =  k1(t)
          y   = [0.04226,0.1593,0.1460,0.1489];
          idx = logical(histc(t,[0,91.25,182.5,273.75,366]));
          y   = y(idx);
      end
      % Nested k2
      function y  =  k2(t)
          y   = [0.008460,0.04959,0.03721,0.04750];
          idx = logical(histc(t,[0,91.25,182.5,273.75,366]));
          y   = y(idx);
      end
  end

Remarks:

Histc uses the following boundary behavior low <= t < up and this is consistent with k(0). In your code k(0) didn't return any value, thus the error you got. In this case it returns 0.04226 (if k1) or 0.008460 (if k2). Note that to include 365 I shifted it to 366 (assuming you'll never have values bigger equal than 366).
Now I get an error on d(1,2) because I don't have code for K, k, gamma etc..
No need to use globals, you can use nested functions which can see parent's workspace.

MATLAB: Changing a parameter in a differential equation

Hi Sean,

you can pass the value of epsilon into the ode function with

[tSol,YSol]=ode45(@(t,Y) quantumsystem(t,Y,epsilon),tRange,Y0);

and

function dYdt = quantumsystem(t,Y,epsilon);
    x=Y(1);
    y=Y(2);
    z=Y(3);
    gamma=1/2;
    c=1;
    dtheta=gamma*c/(t^2 + c^2);
    dxdt=(-1/epsilon)*y - dtheta*z; 
    dydt=(1/epsilon)*x;
    dzdt=dtheta*x;
    dYdt=[dxdt;dydt;dzdt];
end

Since ode45 picks the time increments, different values of epsilon may lead to different lengths of tSol and YSol, but the information is there. One of the easier ways to deal with that issue is to set up a for loop with different values of epsilion and use

plot(tSol,z)
hold on

inside the loop as you are doing, only without the color specification.

Best Answer

Related Solutions

MATLAB: Variable coefficients-help needed please!!!

MATLAB: Changing a parameter in a differential equation

Related Question