Azzera filtri
Azzera filtri

Solving ODE with piecewise driving force

4 visualizzazioni (ultimi 30 giorni)
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?

Risposta accettata

Sam Chak
Sam Chak il 5 Apr 2023
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
  1 Commento
Moosejr
Moosejr il 13 Apr 2023
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')

Accedi per commentare.

Più risposte (0)

Categorie

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

Prodotti


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by