MATLAB Answers

How to determine convergence of ode45

조회 수: 16(최근 30일)
Kostas
Kostas 16 Sep 2019
Commented: Kostas 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

Aquatris
Aquatris 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.
Kostas
Kostas 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).

로그인 to comment.

답변 수 (0)

이 질문에 답변하려면 로그인을(를) 수행하십시오.


Translated by