steady state criterion

Question

Rose 2011-4-6

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/4858-steady-state-criterion

编辑： Vyas 2016-4-26

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Vyas 2016-4-26

编辑：Vyas 2016-4-26

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 2011-4-6

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/4858-steady-state-criterion#answer_6818

在 MATLAB Online 中打开

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Rose 2011-4-6

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?

Rose 2011-4-7

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.

Rose 2011-4-7

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 2011-4-7

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

请先登录，再进行评论。

Answer 2

Laura Proctor 2011-4-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/4858-steady-state-criterion#answer_6813

在 MATLAB Online 中打开

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.