Obtaining value of highest order derivative over simulation period

I am working with ode45 and my system of equations is a 14 dof second order differential equation system. Now on simulation MATLAB returns the value upto first derivative. I want to obtain the value of second order derivatives, how to go about that?

(In some places people have answered to find second order derivative using the array obtained for variable and their first derivatives but in my case the second order derivative is function of other second other derivatives so that is not possible.)

[Torsten]

If you use ODE45, your system must be solved explicitly for the 2nd order derivatives. So do as suggested: call the function in which you prescribe the derivatives to get the 2nd order derivatives.

If your system does not allow this, please include it here for inspection.

Answer 1

[John D'Errico]

Can you use schemes like ODE45, combined with gradient to compute exact estimates of the high order derivatives? Sadly, no. Higher order derivatives will have tiny errors in them, as computed. This happens for various reasons, due to numerical precision and the approximations used. To get truly deeply enough into it would require a numerical analysis course, probably at a grad school level. But let me try a simple example. As my example, I'll use a simple second order ODE, with only ONE variable. Off the top of my head, lets pick this one:

y'' = y' - 4*y

y(0) = 1

y'(0) = 2

I just picked it out of a hat. A VERY simple second order ODE. First, a symbolic solution to compare things to.

syms y(t)
dy = diff(y);
yexact = dsolve(diff(y,2) == dy - 4*y,y(0) == 1,dy(0) == 2)

yexact = 

Now I'll convert this to a problem that ODE45 can solve. We do so by creating a second variable for the first derivative of y, so we turn it into a system of differential equations, each of which is first order. As such, the problem becomes

y' = z

z' = z - 4*y

y(0) = 1

z(0) = 2

This is perfectly standard. In fact, I prefer doing it by hand to using a tool like odeToVectorField, because if I rely on tools to do ALL of my mental work for me, then I get lazy. Having said that, I have limited time and space to do this here, so I'll use odeToVectorField anyway, since that way I'll get it right the first time. My god, I am such a hypocrite. ;-)

[V,Y] = odeToVectorField(diff(y,2) == dy - 4*y)

V = 

Y = 

M = matlabFunction(V,'vars', {'t', 'Y'})

M = function_handle with value:
    @(t,Y)[Y(2);Y(1).*-4.0+Y(2)]

tspan = [0,10];
y0 = [1 2];
[t45,y45] = ode45(M,tspan,y0);

plot(t45,y45(:,1),'bo')
hold on
fplot(yexact,[0,10],'r-')
hold off

In fact, the ODE45 solution should overlay nicely on top of the exact solution. You might think it is exact, but it is not of course, but only an approximation itself. We can see that by looking at the differences.

yexactfun = matlabFunction(yexact)

yexactfun = function_handle with value:
    @(t)exp(t./2.0).*cos((sqrt(1.5e+1).*t)./2.0)+(sqrt(1.5e+1).*exp(t./2.0).*sin((sqrt(1.5e+1).*t)./2.0))./5.0

plot(t45, y45(:,1) - yexactfun(t45))

As you can see, there are in fact small differences, ones that grow in size as t grows larger. They don't really show up well on the scale of the original plot itself. You can see the scale of these differences as being on the order of 1%, so small, but not remotely perfectly exact. Remember, that ODE45 is itself an approximative tool. It has tolerances all of its own, and I used the default tolerances in ODE45.

Now, we can extract the second derivative. At least we can try to do so. I can think of several ways to do so. We can use gradient on the first derivative as produced by ODE45. or we can use gradient to approximate the second derivative from y as predicted from ODE45.

y2app1 = gradient(y45(:,2),t45); % apply gradient to y'
plot(t45,y2app1)

y2app2 = gradient(gradient(y45(:,1),t45),t45); % compute the second derivative of y from y, using consecurtitive calls to gradient.
plot(t45,y2app2)

However, you can see both of those gradient calls both have some issues. They are not as smooth as I want them to be. Remember that gradient uses finite differences to compute a derivative. And since what we applied gradient to was itself only an approximation, we have some serious problems. I can get into some depth here, where I explain that the computation of a derivative is something called an ill-posed problem, which amplifies any tiny noises in the signal. If you think of it, look at the formula for a derivative:

dy/dt = (y(t + dt) - y(t))/dt

Now in theory, we take the limit as dt goes to zero, and out drops a derivative. But gradient does not take any limits. Instead, gradient can use only finite difference approximations.

When I called differentiation an ill-posed problem, now I want to talk about Fredholm integral equations of the first kind, probably worth an entire graduate level course itself. Sigh. Can you just take my word for it, that numerical differentiation is a difficult thing to do accurately?

The true second derivative is given here:

dy2fun = matlabFunction(diff(yexact,2));
y2exact = dy2fun(t45);
plot(t45,y2exact)

So the use of gradient to compute the second derivative is a bad one, because it itself uses a finite difference approximation, and the result of ODE45 is computed as an approximation, effectively a variety of Runge-Kutta.

A better way to compute the second derivative is you just use your function itself, thus M. Unfortunately, M is not itself fully vectorized, as produced by matlabFunction. So I'll rewrite it in vectorized form.

M

M = function_handle with value:
    @(t,Y)[Y(2);Y(1).*-4.0+Y(2)]

Mvect = @(t,y) [y(:,2),-4*y(:,1)+y(:,2)];
dy45 = Mvect(t45,y45);
plot(t45,dy45(:,2))

And that is clearly a decent prediction of the second derivative, as opposed to that which gradient produces. Will it be exact? NOOOOOO!!!!!!!

plot(t45,dy45(:,2) - y2exact)

Again, even these second derivatives predictions are not exactly correct. Remember that ODE45 is an APPROXIMATION itself.

This rabbit hole can get very deep, Can I stop now? Please?

Answer 2

[Star Strider]

One option would be to use the gradient function on the second column of the differential equation output (the integrated function value) as many times as necessary to get the highest derivative you want.

dydt = gradient(y(:,2)) ./ gradient(t);

d2ydt2 = gradient(dydt) ./ gradient(t);

(assuming independent variable ‘t’) and so for the rest.

[Shantanu Chhaparia]

Hi, thanks for the answer. I tried using gradient function, and have a followup question. When I compared the first derivative array that ode45 returned and the first derivative that I calculated by using gradient(yValues)./gradient(tValues) they are coming to be slightly different. Although the difference is two or three order less than the values of the derivative but shouldnt it be exactly same?

And is there some way in which we can make changes to the ode45 call, so that it returns the second derivative as well. Actually I have a mechanical system and need to further calculate the interaction forces. (which can be done by separately putting the array obtained but asking if there is any concise way other than that)

Thanks

[John D'Errico]

NO! Gradient will not produce EXACTLY the same prediction as that what you will get from other methods. That is because gradient produces an APPROXIMATION itself, based on a simple finite difference scheme. So it cannot give you exactly the same result. In fact, the higher order derivatives themselves will have some errors in them, because they were computed using a scheme that has approximations in it.

[Star Strider]

@Shantanu Chhaparia — You can get the derivatives as your function calculates them directly from the function. It requires a loop and the solved system:

[t,y] = ode45(odefcn, tspan, iinitconds);

for k = 1:numel(t)

derivs(:,k) = odefcn(t(k),y(k,:));

end

.

[John D'Errico]

Again, as my answer shows, even this is NOT exact!!!! It will be far better than the use of gradient, which is a sadly terrible idea, even if it seems like it could be ok.

[Star Strider]

It does not have to be exact, just close enough.

[John D'Errico]

My point was that @Shantanu Chhaparia observed differences, and did not understand what was happening. And while you may decide it needs only to be close enough, that is arguably not your decision to make, as to whether any differences are small enough.