How to output derivative from ODE 45

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.

Answer 1

[Reen]

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.

[stuckedMD]

Thank you. I can now return the derivatives. But, why is y(t) doesn't equal with y(t-1) - dy(t-1)?

(Apologies if this is a basic math problem!)

[Walter Roberson]

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]

@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.

Answer 2

[Star Strider]

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

[stuckedMD]

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]

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.