How to output derivative from ODE 45

49 次查看（过去 30 天）
stuckedMD
stuckedMD 2017-9-7
评论： James Tursa 2017-9-7
I have the following function:
function dy = fluboard(t, pop, b, d, N)
dpop(1) = -b*pop(1)*pop(2)/N;
dpop(2) = b*pop(1)*pop(2)/N - d*pop(2);
dpop(3) = d*pop(2);
dy = [dpop(1); dpop(2); dpop(3)];
And call it by:
h=0.01; tspan=[0 20]; pop0=[760 3 0]; b=1.66; d=1/2.2; N=763;
[t y] = ode45(@fluboard, tspan, pop0, h, b, d, N);
How can I get the derivatives at each time point? Thanks.

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

回答（2 个）

Reen
Reen 2017-9-7
If you want other varialbes along with the derivative, you can modify the function you call using ode45 to:
function [dy other_variable]= fluboard(t, pop, b, d, N)
If you just want the derivative you can keep it the same. Now you can run through the code you get your y vector. After that, run through a for loop calling the fluboard function to get your desired output. It should look something like this:
for i=1:length(y)
[dy(i) other_variable(i)] = fluboard(t(i), pop(:,i), b, d, N); % may have to use pop(i,:) instead depending on how the matrix is set up
end
That will just run through your function to find the derivative at each point.
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2017-9-7
The ode routines do not evaluate the function at each time point: they evaluate at nearby time points and predict the value at the specific time points. Any given output might or might not have been evaluated exactly.
James Tursa
James Tursa 2017-9-7
@stuckedMD: If y is linear, then you would get y(t) = y(t-1) + dy(t-1). But y is not linear, so this relationship does not hold. That's the whole reason for calling the ode routines in the first place, because you are trying to solve problems where the solution is not simple and linear.

请先登录，再进行评论。

Star Strider
Star Strider 2017-9-7
编辑：Star Strider 2017-9-7
Taking the numerical derivatives using the gradient function is the easiest solution. For best results, this requires that tspan be a vector with a constant sampling interval, so I added that.
My approach:
fluboard =@(t, pop, b, d, N) [-b*pop(1)*pop(2)/N; b*pop(1)*pop(2)/N - d*pop(2); d*pop(2)];
h=0.01; tspan=[0 20]; pop0=[760 3 0]; b=1.66; d=1/2.2; N=763;
tspan = linspace(0,20);
[t y] = ode45(@(t,pop)fluboard(t, pop, b, d, N), tspan, pop0);
[~,dy] = gradient(y, mean(diff(tspan)));
figure(1)
subplot(2,1,1)
plot(t, y)
title('Solution')
grid
subplot(2,1,2)
plot(t, dy)
title('Derivatives')
grid
  2 个评论
stuckedMD
stuckedMD 2017-9-7
Your suggestion also worked. Thanks. But I'm confused about the relationships between y and dy. I thought y(t) should be equal with y(t-1) - dy(t-1), but it is not. What am I missing? Also, in a different occasion, using constant time interval caused an error. Can this be prevented?
Star Strider
Star Strider 2017-9-7
My pleasure
The gradient function calculates the central difference numerical derivative, except at the edges or ends, where it calculates a one-sided derivative. See the documentation on gradient for a full description.
You did not describe the error, so I cannot describe a way to prevent it. My code as I posted in my Answer ran without error.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by