Using a another function within ODE45

23 visualizzazioni (ultimi 30 giorni)
I have a function that is part of an ODE, however I have no idea how to get it to work with the ODE45 function. I basically want to pass m_t into disp_vel (at the bottom), for values of t. I have also tried using an anonymous function, however this doesn't seem to work either.
The function m_f should be linear, with grad = fluid_r, until it reaches m_f_max, at which point it is constant.
function HarmonicMotion
%Default Values
m_c = 2; % Container Mass
s1 = 16; % Spring Constant 1
s2 = 20; % Spring Constant 2
c = 2; % Damping Coef
d = 0.6; % Gap
t_total = 10; % Total Integration Time
dt = 0.001; % Integration step
xv0 = [2 0]; % Init Displacement & velocity
tvals = 0:dt:t_total; % A 1D array for the time values
m_f_max = 4; % Max mass of fluid
fluid_r = 1; % Fluid mass rate
data = [m_c, s1, s2, c, d, m_c, fluid_r]; % Setting Array for input into functions
function plot_wf(data) % plots a function with the effects of fluid included
[t xv2] = ode45(@(t,x) calc_wf(x, data), tvals, xv0);
figure
plot(t, xv2(:,:),[0 t_total], [d d], 'k-.');
xlabel('time');legend('Displacement / m', 'Velocity / ms^{-1}', 'Gap');grid;
function [disp_vel] = calc_wf(xv, data) % Calcs values for the case a fluid;
x = xv(1); v = xv(2); %1st Col = x, 2nd = vel;
i = 1;
% m_f(i) = 0;
% for t = 0:10*dt:4;
% if m_f <= m_f_max;
% m_f(i+1) = m_f(i) + (fluid_r*(10*dt));
% end
% i = i + 1;
% end
if x > d %Same as for calc_nf
s2_x = x +d;
else
s2_x = 0;
end
disp_vel = [v; -(s1*x + c*v + s2*s2_x)/m_t(1,1)]; % HOW TO GET THIS TO RESPOND TO DIFFERENT VALUES OF m_t
end
end
end
  2 Commenti
Walter Roberson
Walter Roberson il 14 Gen 2021
m_f = @(t)
Invalid syntax. You need to complete the anonymous function
m_t = m_c + m_f
m_f is a function handle and must be invoked to give a result.
Fawaz Zaman
Fawaz Zaman il 14 Gen 2021
Sorry, I should've edited that part out (I was just testing it out), even when it was completed, it would still not work with ODE45. I've corrected this in the question, thank you.

Accedi per commentare.

Risposta accettata

Alan Stevens
Alan Stevens il 14 Gen 2021
More like this (but note the comments near the end):
%Default Values
m_c = 2; % Container Mass
s1 = 16; % Spring Constant 1
s2 = 20; % Spring Constant 2
c = 2; % Damping Coef
d = 0.6; % Gap
t_total = 10; % Total Integration Time
dt = 0.001; % Integration step
xv0 = [2 0]; % Init Displacement & velocity
tvals = 0:dt:t_total; % A 1D array for the time values
m_f_max = 4; % Max mass of fluid
fluid_r = 1; % Fluid mass rate
data = [m_c, s1, s2, c, d, fluid_r ,m_f_max]; % Setting Array for input into functions
[t, xv2] = ode45(@(t,x) calc_wf(t, x, data), tvals, xv0);
figure
plot(t, xv2(:,:),[0 t_total], [d d], 'k-.'), grid
xlabel('time');
legend('Displacement / m', 'Velocity / ms^{-1}', 'Gap')
figure
plot(t,min(fluid_r*t + m_c, m_f_max)), grid
xlabel('time'),ylabel('fluid mass')
function [disp_vel] = calc_wf(t, xv, data) % Calcs values for the case a fluid;
% Extract data
m_c = data(1); s1 = data(2); s2 = data(3); c = data(4);
d = data(5); fluid_r = data(6); m_f_max = data(7);
% Extract displacement and velocity
x = xv(1); v = xv(2); %1st Col = x, 2nd = vel;
% Mass at time t
m_t = min(fluid_r*t + m_c, m_f_max);
if x > d %Same as for calc_nf %%%%%%%%%% In calc_nf you had x<-d
s2_x = x + d; %%%%%%%%%%%%% Do you mean this?
%%%%%%%%%%%%% Ifx and d are both positive
%%%%%%%%%%%%% shouldn't you have x-d here?
else
s2_x = 0;
end
disp_vel = [v; -(s1*x + c*v + s2*s2_x)/m_t];
end
  5 Commenti
Walter Roberson
Walter Roberson il 14 Gen 2021
Global is recommended against, as it is the slowest form of variable and can be the most difficult to debug.
With respect to event functions, I recommend the ballode example.
Alan Stevens
Alan Stevens il 14 Gen 2021
1. Try
doc ode event location
2. Best to avoid global variables in general. See https://matlab.fandom.com/wiki/FAQ#Are_global_variables_bad.3F for example.

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by