Solving ODE with piecewise driving force

I have the following equation of motion

with the driving force

I want to solve the equation numerically, as I might want to slightly alter the driving term later.

I can solve the ODE with by using ODE45. However, If I now try to include the condition that when the driving force reaches zero it should stay zero, I get the following:

tspan = [0 5];
V0 = 0;

[t,V] = ode45(@(t, V) Force(t,V), tspan, V0);
f = Force(t,V);

figure(1)
subplot(1,2,1)
plot(t,V)
xlabel('t')
ylabel('y')

subplot(1,2,2)
plot(t,f)
xlabel('t')
ylabel('Driving Force')

function RHS = Force(t,V)
RHS = 2*exp(-t) - V;
if RHS < 0
    RHS = 0;
end
end

The solution y vs t looks OK, in the sense that the object stops being accelerated when the driving force reaches zero. However, given what I have written in the force function I would expect the driving force to become zero. The plot shows that this is clearly not the case, which also makes me suspicous of the solution for y vs t.

My question is then, how do I correctly solve the ODE with the piecewise driving force term?

Answer 1

[Sam Chak]

Hi @Moosejr

In this case, you can use the deval() function. I also fixed the conditional statement in the ODE.

tspan   = [0 5];
V0      = 0;
sol     = ode45(@(t, V) Force(t,V), tspan, V0);
t       = linspace(0, 5, 5001);
[V, Vp] = deval(sol, t);   % V is output, Vp is V-prime (V')

figure(1)
subplot(2,1,1)
plot(t, V, 'linewidth', 1.5), grid on
xlabel('t')
ylabel('y')

subplot(2,1,2)
plot(t, Vp, 'linewidth', 1.5), grid on
xlabel('t')
ylabel('Driving Force')

function RHS = Force(t, V)
    fcn = 2*exp(-t) - V;
    if fcn > 0
        RHS = 2*exp(-t) - V;
    else
        RHS = 0;
    end
end

[Moosejr]

Thank you.

For future reference to plot the RHS of the ODE directly instead of the first derivative one can write:

f = [];
for i = 1:length(t)
    f = [f,Force(t(i),V(i))];
end
subplot(3,1,3)
plot(t,f)
xlabel('t')
ylabel('Driving Force')