steady state criterion

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.

[Vyas]

Hello,

I do have a similar set of ODE's which solves a chemical kinetics system. I am using some constraints to solve these equations and I am able to solve for few set of unknowns. However when I use more unknowns (~8) with 2 constraints to solve for rest 6 unknowns using ODE, I see the problem will not be solved till the final time. For e.g., I have specified initial time and final time as [t0 tf] which i pass to the ODE15s solver. So for the higher number of unknowns, the ODE solver will not solve until 'tf' has achieved (not blowing up as well). Can anyone help me why ODE solves only to a fraction of final time 'tf'? And also tell me how to fix this? I have seen events being used to stop the code when steady state achieved. I want something reverse. I want to know why the ODE solver do not finishes till the specified time!!

Answer 1

[Matt Tearle]

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);

[Andrew Newell]

You'd need to check that the norm is not less than the tolerance before the integration starts.

[Rose]

Thanks a lot! I was so stressed but now I think it is possible to solve this problem. Thanks once again! I like your idea and I tried to implement it but it didn't work. I am posting my whole code. May be you can help me out!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This is my subroutine which I am using to state my ODEs:

function dy = funl1(t,y)

global k1 k2 k3 mu d1 d2 r K k alpha gamma n;

dy(1,1) = (k1.*(y(4)-y(1)).*y(3))./(y(2)+y(3)-1e-12) - k2.*y(1).*(y(2)./(y(2)+y(3)-1e-12)) ;

dy(2,1) = (mu.*(y(2).^n)/(K^n+y(2).^n)).*exp((-y(1)).*k)-(k3.*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d1+gamma.*y(4)).*y(2);

dy(3,1) = (k3.*y(1).*y(2))./(y(2)+y(3)-1e-12)-(d2+gamma.*y(4)).*y(3);

dy(4,1) = r.*y(4).*(1-(y(4)./(alpha.*(y(2)+y(3)-1e-12))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is the main code :

global k1 k2 k3 mu d1 d2 r K k alpha gamma n;

ks =[0.04226 0.008460 0.03384 0.00 0.005263 0.005300 0.0045 10000 0.00003125 0.4784 10^(-6) 2];

k1 = ks(1);

k2 = ks(2);

k3 = ks(3);

d1 = ks(5);

d2 = ks(6);

mu = ks(4);

r = ks(7);

K = ks(8);

k = ks(9);

alpha = ks(10);

gamma = ks(11);

n = ks(12);

finaltime = 10000;

tspan = [0 finaltime];

y0=[0;10000;0;0];%winter

options = odeset('RelTol',1e-12,'AbsTol',1e-12);

sol = ode15s(@funl1,tspan,y0,options);

X = sol.x;

Y = (sol.y)';

plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3), X, Y(:,4))

xlabel('time')

ylabel('Y')

legend('m', 'X', 'y', 'M')

figure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I don't know what this copy and paste thing? Do I need to write my equations again?

[Matt Tearle]

Yes, just copy and paste from funl1 to your event function. Or, for a more elegant solution, call your ode function (assuming scope/visibility allows it).

While I'm here, let me also say: don't use global! Global is the last resort of the wicked. A better way:

%%%%% funl1.m %%%%%

function dy = funl1(~,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n)

dy(1,1) = [etc]

%%%%% odeevent.m %%%%%

function [x,isterm,dir] = odeevent(~,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n)

dy = funl1([],y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n);

x = norm(dy) - 1e-5;

isterm = 1;

dir = 0;

%%%%% main script %%%%%

ks =[0.04226 0.008460 0.03384 0.00 0.005263 0.005300 0.0045 10000 0.00003125 0.4784 10^(-6) 2];

k1 = ks(1);

k2 = ks(2);

[etc]

odefun = @(t,y) funl1(t,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n);

efun = @(t,y) odeevent(t,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n);

options = odeset('Events',efun);

sol = ode15s(odefun,tspan,y0,options);

odefun and efun are anonymous function handles; they are essentially wrapper functions to funl1 and odeevent with all the parameters embedded in them (just t and y left variable).

[Rose]

I tried this thing, but it is giving me error- input argument k1 is undefined.

%%%%%%%funl2.m%%%%%%%%%%%%

function dy = funl2(~,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n)

dy(1,1) = (k1.*(y(4)-y(1)).*y(3))./(y(2)+y(3)+1e-12) - (k2*y(1).*y(2))./(y(2)+y(3)+1e-12) ;

dy(2,1) = (mu.*(y(2).^n)./(K.^n+y(2).^n)).*exp(-y(1).*k)-(k3*y(1).*y(2))./(y(2)+y(3)+1e-12)-(d1+gamma.*y(4)).*y(2);

dy(3,1) = (k3*y(1).*y(2))./(y(2)+y(3)+1e-12)-(d2+gamma.*y(4)).*y(3);

dy(4,1) = r.*y(4).*(1-(y(4)./(alpha.*(y(2)+y(3)+1e-12))));

%%%%%%%%%%%%%odeevent.m%%%%%%%%%%%%%%%

function [x,isterm,dir] = odeevent(~,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n)

dy = funl2([],y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n);

x = norm(dy) - 1e-5;

isterm = 1;

dir = 0;

%%%%%%%%%%%steady_state.m%%%%%%%%%%%

ks =[0.04226 0.008460 0.03384 0.00 0.005263 0.005300 0.0045 10000 0.00003125 0.4784 10^(-6) 2];

k1 = ks(1);

k2 = ks(2);

k3 = ks(3);

mu = ks(4);

d1 = ks(5);

d2 = ks(6);

r = ks(7);

K = ks(8);

k = ks(9);

alpha=ks(10);

gamma=ks(11);

n = ks(12);

finaltime = 2900;

tspan = [0 finaltime];

options=odeset('RelTol',1e-6,'AbsTol',...

1e-12,'Stats','on');

y0 = [0;20000;0;0];

odefun = @(t,y) funl2(t,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n);

efun = @(t,y) odeevent(t,y,k1,k2,k3,mu,d1,d2,r,K,k,alpha,gamma,n);

options = odeset('Events',efun);

sol = ode15s(@funl2,tspan,y0,options);

X = sol.x;

Y = (sol.y)';

plot(X,Y(:,1),'r-', X,Y(:,2),'g-', X,Y(:,3),'b--', X,Y(:,4))

xlabel('time')

ylabel('Y')

legend('m', 'X', 'y', 'M')

figure

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Can you please check this. I don't know why is not working with me.

[Matt Tearle]

You're still using @funl2 as the input to ode15s, instead of the new wrapper function odefun

sol = ode15s(odefun,tspan,y0,options);

[Rose]

now the error is:

too many input arguments!!!

[Rose]

error using ====> funl2

[Matt Tearle]

I don't get any error. I copy/pasted exactly what you have above, and the only change I made was to replace "@funl2" with "odefun". It works fine.

I suggest: using the error message to trace back; using "which" to make sure you're getting the files you expect; making sure your functions are all saved (ie no unsaved changes); using the debugger.

The error message suggests a mismatch between the argument lists in the function definition and where you call it. What you have above works fine for me, so perhaps something was changed or you have more than one version of funl2 (eg an old one) and it's finding that instead of the one you want.

[Rose]

It worked!!!! Thanks Thanks Thanks Matt :)

But I think it is not taking into account the RelTol and AbsTol

How can I fix it?

[Rose]

I think I did it :) Thanks a lot for your help! I am so relaxed now!!

Answer 2

[Laura Proctor]

You can set the error tolerance of the ODE solver using the function ODESET and then pass in the structure created from ODESET into your ODE solver.

For example:

options = odeset('RelTol',1e-2,'AbsTol',[1e-3,1e-3,1e-4,1e-2]);
[T,Y] = ode15s(@myFcn,t,y0,options);

Try this out to see if you're able to get better results.

You may also be interested in event handling, which is where you can make your simulation stop once an event happens (such as when the output goes below a specific value). You also do this using ODESET function - look at the topic Event Location Property in the documentation for ODESET for more details on this.

[Rose]

I tried it and it gave me better results but I want to apply the steady state criterion particularly!!