How to determine convergence of ode45

17 views (last 30 days)
KostasK
KostasK on 16 Sep 2019
Commented: KostasK on 18 Sep 2019
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 Comments
Aquatris
Aquatris on 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.
KostasK
KostasK on 18 Sep 2019
Yes, this specific problem is simple enough such that this approach could suffice. I think for my above question, its better to calculate the Residuals than the error convergence (which should happen anyway since the solver does not return any warnings/errors).

Sign in to comment.

Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!