Using a another function within ODE45

24 Ansichten (letzte 30 Tage)
Fawaz Zaman
Fawaz Zaman am 14 Jan. 2021
Kommentiert: Alan Stevens am 14 Jan. 2021
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 Kommentare
Walter Roberson
Walter Roberson am 14 Jan. 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 am 14 Jan. 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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Alan Stevens
Alan Stevens am 14 Jan. 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 Kommentare
Walter Roberson
Walter Roberson am 14 Jan. 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 am 14 Jan. 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by