Obtaining value of highest order derivative over simulation period

15 views (last 30 days)
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.)
  1 Comment
Torsten on 16 May 2022
Edited: Torsten on 16 May 2022
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.

Sign in to comment.

Answers (2)

John D'Errico
John D'Errico on 17 May 2022
Edited: John D'Errico on 17 May 2022
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:
tspan = [0,10];
y0 = [1 2];
[t45,y45] = ode45(M,tspan,y0);
hold on
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:
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'
y2app2 = gradient(gradient(y45(:,1),t45),t45); % compute the second derivative of y from y, using consecurtitive calls to gradient.
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);
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 = function_handle with value:
Mvect = @(t,y) [y(:,2),-4*y(:,1)+y(:,2)];
dy45 = Mvect(t45,y45);
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?

Star Strider
Star Strider on 16 May 2022
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.
John D'Errico
John D'Errico on 17 May 2022
Edited: John D'Errico on 17 May 2022
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.

Sign in to comment.


Community Treasure Hunt

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

Start Hunting!