ODE45 stop when steady state occurs in a periodic function

Question

KostasK 2019-12-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/497192-ode45-stop-when-steady-state-occurs-in-a-periodic-function

回答： Vitaly Kheyfets 2022-9-22

Hi all,

I have a typical set of equations for a forced spring-mass-damper system, which I have managed to solve successfully. My problem is that I would like to know how is it possible to stop the equations when steady state is reached in a periodic response.

I know in non-periodic responses this task is trivial, as I can provide a criterion through an event function that wen achieved, it stops the integration. For example the criterion for some response

can be when

.

However for the case of a periodic function this criterion has to depend on the results of the previous iterations, such that for a periodic response

with a period T, i can set a criterion such as

(where I choose C). Since I don't have access to previous iterations, how would I do this in this case?

Below is an example of the type of response that I am looking at.

Thanks for your help in advance,

KMT.

clear ; clc
% Inputs
    N = 3 ; 
    j = 50 ; 
    k = 200 ; 
    c = 50 ;
    cc= 50 ; 
    f = 100 ; 
    om = 5 ; 
    tmax = 20 ;
    dth0 = 5 ; 
% Matrices
    J = diag(repmat(j, N, 1)) ; 
    K = tridiag(k, N) ; 
    C = tridiag(c, N) ; 
    CC= diag(repmat(cc,N, 1)) ; 
    F = diag(repmat(f, N, 1)) ; 
    P = linspace(0, 2*pi*(1 - 1/N), N)' ; 
% Solution
    tspan = [0 tmax] ; 
    th0 = [zeros(N, 1) ; repmat(dth0, N, 1)] ; 
    options = odeset('Events', @(t, th) eventfcn(t, th, N), 'RelTol', 1e-4) ;
    
    [t, Mth, te, Mthe] = ode45(@(t, th) odecfn(t, th, J, K, C, CC, F, P, N, om), tspan, th0, options) ;
% Plotting
    figure
    plot(Mth(:,N), Mth(:,N+1:2*N))
    hold on
    plot(te, Mthe(:,N+1:2*N), 'o')
    grid on
% Events Function
function [va, is, di] = eventfcn(~, th, N)
va = th(N) < 5 ;
is = 0 ; 
di = 0 ;
end
% ODE Function
function dthdt = odecfn(t, th, J, K, C, CC, F, P, N, om) 
% Indexes
i = N ; 
j = N+1 ; 
k = 2*N ; 
% Preallocate
dthdt = zeros(k, 1) ; 
% Equation
dthdt(1:i) = th(j:k) ; 
dthdt(j:k) = J \ (-[C+CC] * th(j:k) - K * th(1:i) + F * (1 + sin(om*t - P))) ; 
end
% Tridiagonal Matrix Function
function M = tridiag(m, N)
M1 = 2*diag(repmat(m, N, 1)) ;
M2 = diag(repmat(m, N-1, 1), -1) ;
M3 = diag(repmat(m, N-1, 1), 1) ;
M = M1 - M2 - M3 ;
M(1,1) = M(1,1) - m ;
M(end,end) = M(end,end) - m ;
end

5 个评论
显示 3更早的评论隐藏 3更早的评论

J. Alex Lee 2019-12-18

在 MATLAB Online 中打开

Oops, I was on a different computer, I was actually running 2017b, which will return error on

sign(true)

Your eventfcn is returning a logical (the result of th(N)<5), but it expects a signed double; it looks for the expression crossing 0, and it looks for conditions from left or right, so it tries to take the sign. I gues TMW reconsidered that it should let sign(true) return 1. In any case, you probably want to fix that line for expected behavior.

So i removed the event and I see your periodicity. I guess you want to stop around t=20?

Regarding the actual question, no I don't mean having 2 ODE functions, something like

tspan = [t0,t1];
% th0 from above
sol = ode45(@(t, th) odecfn(t, th, J, K, C, CC, F, P, N, om), tspan, th0, options) ;
% do whatever you want with sol to check for periodicity
tspan = [sol.x(end),sol.x(end) + (t1-t0)];
th0 = sol.y(end,:); % or (:,end), I forget which
sol_b = ode45(@(t, th) odecfn(t, th, J, K, C, CC, F, P, N, om), tspan, th0, options) ;

One problem is how to combine sol and sol_b to interpolate on the whole solution...it might be as simple as concatenating the fields as below, but i'm not sure exactly how deval() works

sol_full.x = cat(sol.x,sol_b.x);
sol_full.y = cat(sol.y,sol_b.y); % again can't remember shape of y, or may depend on how you set it up
sol_full.solver = sol.solver; % so you'd have to make sure you use ode45 for both

As for how to actually check for periodicity once you have access to a solution history, do you have that sorted out, because I'm curious to know how you're doing that part.

Walter Roberson 2019-12-18

Note that boolean conditions never cross 0, as they can only assume the value 0 (false) or 1 (true) and never negative.

KostasK 2019-12-18

编辑：KostasK 2019-12-18

I see what you mean. I was hoping that I was going to be able to do something like obtain successive periods of my response, and determine the relative errors between those to see how small it gets and based on that determine weather I have reached 'steady state'. This is of course something that I can do after I have access to my solution, having had ODE45 do its job first, but it is increasingly looking like its not something that I can do whilst the solver is running.

Thanks for the help :)

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2019-12-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/497192-ode45-stop-when-steady-state-occurs-in-a-periodic-function#answer_407006

Any attempt to do computations on differential equations in which the results depend upon the computations at a different earlier time, must be coded as Delay Differential Equations, probably using dde23() . ddeset() can be used to add Event functions, which you could use to signal termination.

Be careful not to test just a single delay, (t) vs (t-period) : during a phase of oscillations it would not be uncommon to find some time at which f(t) == f(t-period) . But perhaps testing a couple of derivatives would be enough. Or code in two lags, so that you have available data for (t), (t-period), (t-2*period)

4 个评论
显示 2更早的评论隐藏 2更早的评论

Walter Roberson 2019-12-20

I guess for this question, it has to be true that DDEs always result in periodic solutions and that the set of lags is somehow exactly related to the period of the solutions...is this the case?

In the general case, DDE are seldom periodic. Consider for example a simple situation in which an input force of water has to flow over a section of ground before falling into something to provide motive power: the output power is delayed relative to the input rather than there being an instant response.

In the original poster's case, we are given that that the steady state of the system being considered is periodic with a known delay. That being the case, we can compare all of the state variables for time (t) to the state variables at time (t-lag): if they differ more than nominally then we are not yet in that steady state. But do not just compare one of the state variables such as output position, since any given position can be the same between two times even though in one case the object is moving up and in the other case the object is moving down. Either compare all of the state variables or else compare fewer state variables at multiple periods.

Note that comparing state variables to the same variables from a previous period will not necessarily work well with chaotic or quasi-periodic systems.

KostasK 2019-12-20

Hi Walter, in regards to your first point, its a yes: in this question the DDEs always result in periodic solution. To be more precise as you said they don't exactly result in periodic solution, but for this case is extremely close.

Walter Roberson 2019-12-20

mass-spring systems are notorious for being sensitive to external vibration that can lead to some fairly notable movement even when the system starts from "rest". But that depends upon the arrangements and the amount of friction in the system. I gather there is established mathematics for determining whether small vibrations will get magnified arbitrarily (it would not surprise me if the calculation was along the lines that eigenvalues with absolute value greater than 1 implied instability relative to small vibrations.)

请先登录，再进行评论。

Answer 2

Vitaly Kheyfets 2022-9-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/497192-ode45-stop-when-steady-state-occurs-in-a-periodic-function#answer_1058985

Hi KMT,

Another option, if you wanted to continue using ode45 and avoid the delay DE, is to simply put the ode45 optoin into a loop. You could run a single cycle in a loop (or 3 cycles) and then check that the solution is periodic. If it did, great, you're done. If it didn't, then you use the final solution of the last iteration as the initial condition of the next loop iteration.

Hope that helps!

Vitaly

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

ODE45 stop when steady state occurs in a periodic function

5 个评论
显示 3更早的评论隐藏 3更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

ODE45 stop when steady state occurs in a periodic function

5 个评论 显示 3更早的评论隐藏 3更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论