ODE45 stop when steady state occurs in a periodic function

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

[J. Alex Lee]

your code errors out (2019b)

Undefined function 'sign' for input arguments of type 'logical'.

Error in odezero (line 46)
    indzc = find((sign(vL) ~= sign(vR)) & (direction .* (vR - vL) >= 0));

Error in ode45 (line 353)
        odezero(@ntrp45,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);

Error in untitled (line 25)
    [t, Mth, te, Mthe] = ode45(@(t, th) odecfn(t, th, J, K, C, CC, F, P, N, om), tspan, th0, options) ;

For clarification, I guess you are not looking for a "steady state", but rather when you have completed some number of cycles of a periodic response.

One change I would make to your code is to return the solution structure from ode45 (use single output argument) so that you have access to the full solution interpolation (using deval).

As for access to solution history: can you just choose an interval, check if the solution so far is periodic to your satisfaction, and if not, re-start the integration from the last point and repeat?

[KostasK]

Hi Alex, thanks for your response.

That's weird you are getting those errors. I just run the code right now on MATLAB Version: 9.7.0.1216025 (R2019b) Update 1, and it runs without errors. Is it possible that I might have some packages that you dont?

On the rest of your response, yes you are correct in that I am looking to see if the periodic response from one cycle matches the next (hence some form of a steady state).

For accesing the solution history by stopping and restarting, do you mean having 2 ode functions, or is it possible to stop and re-start with 1 function?

[J. Alex Lee]

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]

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

[KostasK]

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]

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)

[J. Alex Lee]

I am curious now how to represent this as a DDE (the concept is new to me)...I see on cursory learning that some DDEs can be represented as ODEs, so if the system above belongs to that class, the answer might be straightforward. In particular, 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?

And in Walter's warning to "not test just a single delay", it seems to me that testing a delay in the ddeXX() context means asserting different lag times...if this is true, wouldn't that mean you are changing the problem itself?

And what if there is no analogous DDE? Can any ODE be represented as a DDE for purposes of "tracking" a solution in time?

[Walter Roberson]

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]

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]

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]

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