ODE45 stop when steady state occurs in a periodic function

26 views (last 30 days)
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 Comments
Walter Roberson
Walter Roberson on 18 Dec 2019
Note that boolean conditions never cross 0, as they can only assume the value 0 (false) or 1 (true) and never negative.
KostasK
KostasK on 18 Dec 2019
Edited: KostasK on 18 Dec 2019
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 :)

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 18 Dec 2019
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 Comments
KostasK
KostasK on 20 Dec 2019
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
Walter Roberson on 20 Dec 2019
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.)

Sign in to comment.

More Answers (1)

Vitaly Kheyfets
Vitaly Kheyfets on 22 Sep 2022
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!