How to determine convergence of ode45
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi all,
I was wondering how it is possible to determine convergence of ode45 for a simple damped mass-spring system. So here is the code for the system:
clear
clc
% Inputs
% Degrees of Freedom
N = 3 ;
% Masses (kg)
m = repmat(100, 1, N) ;
% Spring coefficients (N/m)
k = repmat(15, 1, N+1) ;
% Damping coefficients (Ns/m)
c = repmat(5, 1, N+1) ;
% Maximum force (N)
f0 = repmat(100, 1, N) ;
% Phase angle between force (rad)
phi = linspace(0, pi, N) ;
% Frequency of force (Hz)
freq = 0.2 ;
% Time of simulation (s)
tmax = 100 ;
% Mass initial positions (m)
xi = zeros(1, N) ;
% ODE Inputs
% Degrees of Freedom
N = length(m) ;
% Mass Matrix
M = diag(m) ;
% Stiffness Matrix
k1 = diag(k(1:end-1) + k(2:end)) ;
k2 = diag( -k(2:end-1), 1) ;
k3 = diag( -k(2:end-1), -1) ;
K = k1 + k2 + k3 ;
% Damping Matrix
c1 = diag(c(1:end-1) + c(2:end)) ;
c2 = diag( -c(2:end-1), 1) ;
c3 = diag( -c(2:end-1), -1) ;
C = c1 + c2 + c3 ;
% ODE Initial Conditions
% Time
tspan = [0 tmax] ;
% Displacement and velocity
x0 = [ xi zeros(1, N) ] ;
% ODE Options
options_set = odeset('Mass', [eye(N) zeros(N) ; zeros(N) M]) ;
% Solution
[t, xSol] = ode45(@(t, x) odefcn(t, x, K, C, f0, freq, phi, N), tspan, x0, options_set) ;
% Results
x = xSol(:, 1:N) ; % Displacement (m)
xdot = xSol(:, N+1:end) ; % Velocity (m/s)
% ODE Function
function dxdt = odefcn(t, x, K, C, f0, freq, phi, N)
% Indices
m = N + 1 ;
n = N * 2 ;
% Preallocate
dxdt = zeros(n, 1) ;
% ODE Equation
dxdt(1:N) = x(m:n) ;
dxdt(m:n) = - K * x(1:N) - C * x(m:n) + f0' .* sin(2*pi*freq * t + phi') ;
end
Now, as I understand, in order to determine convergence I will have to use an EventsFcn. However what perplexes me is how to express the convergence formula within the odefcn, as I can't see how I can compare the one step with the next.
Kind Regards,
KMT.
2 Kommentare
Aquatris
am 17 Sep. 2019
Since it is a simple mass-spring-damper system, why don't you just check the pole locations (eigenvalues of A matrix, A =[zeros(n) eye(n); -M\K -D\K]). If all poles are on the left half plane (negative real parts) than the ODE will converge.
Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!