ODE45 stop when steady state occurs in a periodic function
7 ビュー (過去 30 日間)
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
N = 3 ;
j = 50 ;
k = 200 ;
c = 50 ;
cc= 50 ;
f = 100 ;
om = 5 ;
tmax = 20 ;
dth0 = 5 ;
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)' ;
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) ;
plot(te, Mthe(:,N+1:2*N), 'o')
% 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)
i = N ;
j = N+1 ;
k = 2*N ;
dthdt = zeros(k, 1) ;
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 ;
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)
その他の回答 (1 件)
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!