MATLAB Answers

How to extract the value of dydt from ode45 function

조회 수: 17(최근 30일)
Kostas
Kostas 13 Sep 2019
Hi all,
I have written a code of a damped mass-spring system, that takes as inputs the relevant mass, damping, stiffness and excitation force matrices/vectors, and outputs the displacement and the velocity of each mass.
What I would like to do is to be able to extract the acceleration as well (i.e. the value of dxdt(m:n) as written in the odefcn function at the bottom of the code) in addition to the displacement and velocity outputs.
Now I know that I can find the gradient of the velocity, or reconstruct my equation of motion using my solutions for the displacement and velocity and solve for the acceleration. However, having tried, I would like to avoid these approaches, as the actual problem that I am being faced with (which has significantly different equations of motion) has noisy velocity vectors, and the equation is quite complicated to reconstruct and solve for the acceleration (doable but I'd rather avoid it).
So her is the code:
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) ] ;
% Solution
[t, xSol] = ode45(@(t, x) odefcn(t, x, M, K, C, f0, freq, phi, N), tspan, x0) ;
% Results
x = xSol(:, 1:N) ; % Displacement (m)
xdot = xSol(:, N+1:end) ; % Velocity (m/s)
% ODE Function
function dxdt = odefcn(t, x, M, 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) = - M \ K * x(1:N) - M \ C * x(m:n) + M \ f0' .* sin(2*pi*freq * t + phi') ;
end
Thanks for your help in advance,
KMT.

  댓글 수: 0

로그인 to comment.

채택된 답변

James Tursa
James Tursa 13 Sep 2019
James Tursa 님이 편집함. 13 Sep 2019
Why can't you just call your odefcn( ) function with your solution t and xSol vector elements as inputs (e.g., in a loop)? Doesn't this give you the accelerations you want?

  댓글 수: 2

Kostas
Kostas 14 Sep 2019
I didn't know I could use the odefcn like that. That's what I am trying to do now, thanks!
John D'Errico
John D'Errico 14 Sep 2019
Yes. It returns the derivative that you wanted anyway. There is another trick you could use, making the acceleration itself a variable, then you would need to extend the set of ODE variables to solve for. More complicated, certainly, and not worth the effort. It would also slow down the code somewhat, and it would force you to modify your code. As I said, not worth the effort when a simple solution exists.

로그인 to comment.

More Answers (0)

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


Translated by