Using a another function within ODE45

조회 수: 24 (최근 30일)
Fawaz Zaman
Fawaz Zaman 2021년 1월 14일
댓글: Alan Stevens 2021년 1월 14일
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
Walter Roberson
Walter Roberson 2021년 1월 14일
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 2021년 1월 14일
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.

댓글을 달려면 로그인하십시오.

채택된 답변

Alan Stevens
Alan Stevens 2021년 1월 14일
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
Walter Roberson
Walter Roberson 2021년 1월 14일
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 2021년 1월 14일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by