清空筛选条件
清空筛选条件

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 中查找有关 Programming 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by