ODE45 stop when steady state occurs in a periodic function
조회 수: 26 (최근 30일)
이전 댓글 표시
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
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.
채택된 답변
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
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
2022년 9월 22일
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
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!