Conditionals within ODE45

1 view (last 30 days)
Fawaz Zaman
Fawaz Zaman on 13 Jan 2021
Commented: Fawaz Zaman on 14 Jan 2021
Essentially, I am trying to modify the equation passed to ode45 depending on the value of y it creates. I have attempted to do this using a piecewise function, but to no avail, I get an
Error using symengine
Unable to generate code for
piecewise for use in
anonymous functions.
error.
Is there any work around this, or perhaps another method I can use to achieve this effect.
Also attached is an image that may help explain the actual physics of the problem. Thanks
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.01; % Integration step
x_0 = 2.0; % Init Displacement
tvals = [0:dt:t_total];
x1 = nf_s1(c, m_c, s1, s2, d, dt, t_total, x_0);
plot(tvals, x1)
function [disp_vel] = nf_s1(c, m_c, s1, s2, d, dt, t_total, x_0)
tvals = [0:dt:t_total];
syms y(t);
diff_x = y(t)-d; % the adj displacment for the second spring
is_s2 = piecewise(y(t) <= d, 0, y(t) > d, 1) % turns 'on/off' the effect of the bottom spring
eqn = s1*y(t) + c*diff(y) + is_s2*s2*(diff_x) == -m_c * diff(y,2);
[V] = odeToVectorField(eqn)
M = matlabFunction(V, 'vars', {'t', 'Y'})
[tvals, disp_vel] = ode45(M, tvals, [x_0 0]);
  1 Comment
Jan
Jan on 13 Jan 2021
Just a remark: Matlab's ODE integratore are designed to handle smooth functions only. The jump of the foirces, when the spring gets or looses contact can disturb the step size controller such that the integration stops with an error or the result is dominated by rounding errors.
Use an event function to stop the integration and restart it with the using the final value as initial value for the next interval with the changed parameters.

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 13 Jan 2021
You could try the following simplistic approach
%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
data = [m_c, s1, s2, c, d];
t_total = 10; % Total Integration Time
dt = 0.01; % Integration step
xv0 = [2.0 0]; % Init Displacement & velocity
tvals = 0:dt:t_total;
[t, x1] = ode45(@(t,x) nf_s1(t,x,data),tvals, xv0);
figure
plot(t, x1(:,1),[0 t_total],[-d -d],'k--'),grid
ylabel('displacement')
hold on
yyaxis right
plot(t,x1(:,2))
ylabel('velocity')
xlabel('time')
function [disp_vel] = nf_s1(~, xv, data)
m_c = data(1); s1 = data(2); s2 = data(3); c = data(4); d = data(5);
x = xv(1); v = xv(2);
if x<-d
y = x+d;
else
y = 0;
end
disp_vel = [v;
-(s1*x + c*v + s2*y)/m_c];
end
The step change in y doesn't seem large enough to adversely affect ode45 significantly here.
  1 Comment
Fawaz Zaman
Fawaz Zaman on 14 Jan 2021
Hi, thanks for your solution and with neating up the code a little bit. I'm quite new to MatLab, so I do appreciate the help a lot.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!