ODE45 stop when steady state occurs in a periodic function

7 ビュー (過去 30 日間)
KostasK 2019 年 12 月 18 日
回答済み: 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,
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
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 ;
% 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))) ;
% 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 ;
  5 件のコメント
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 :)



Walter Roberson
Walter Roberson 2019 年 12 月 18 日
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 件のコメント
Walter Roberson
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.)


その他の回答 (1 件)

Vitaly Kheyfets
Vitaly Kheyfets 2022 年 9 月 22 日
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!


Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by