Solving ODE with piecewise driving force

2 次查看（过去 30 天）
Moosejr
Moosejr 2023-4-5
评论： Moosejr 2023-4-13
Ran in:
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?

请先登录，再回答此问题。

采纳的回答

Sam Chak
Sam Chak 2023-4-5
Ran in:
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
  1 个评论
Moosejr
Moosejr 2023-4-13
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')

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by