Second time derivative of ODE45 solution looks bumpy

I'm trying to solve a second order differential equation in matlab using the ODE45 solver. I'm quite new to this and would like to plot the solution, as well as the first and second derivatives of the solution. I'm sure there is a better way to achieve this anyway, but I'm taking the numerical derivative of the ODE45 soution twice and plotting it. This worked for a couple differential equations, but for this specific differential equation produces odd results.

Here is my code:

dt = .001;
t=[0:dt:5];
initz=[0; 0];
[t,z]=ode45(@f3, t, initz);
x1 = gradient(z(:,2))./dt;
x2 = gradient(x1)./dt;
figure (3)
plot(t, z(:,2), t, x1, t, x2)

Here is my function:

function dz=f3(t,z)
% z(1) = x1 = dx/dt;   z(2)= x 
dz=[-100*z(1)+cos(5*t-pi/3); z(1)];
end

Here is the function I'm trying to solve:

And here is the output plot showing the issue I'm having.

If theres a better way of doing this, I'm all ears. But I would also like to know what's going wrong here if anyone can figure it out.

Thank you!

Answer 1

[Torsten]

t=[0 5];
initz=[0; 0];
[t,z]=ode15s(@f3, t, initz, odeset('RelTol',1e-12,'AbsTol',1e-12));
for i=1:numel(t)
  dz = f3(t(i),z(i,:));
  x2(i) = dz(1,1);
end
plot(t, z(:,2), t, z(:,1), t, x2)
ylim([-0.05 0.05])

Here is my function:

function dz=f3(t,z)
% z(1) = x1 = dx/dt;   z(2)= x 
dz=[-100*z(1)+cos(5*t-pi/3); z(1)];
end

[Paul Miller]

Thank you, even just replacing my derivative section with your's fixed my plot despite that I still used ODE45, however I still don't understand why my gradient method produced such poor results. If you're able to explain what I was doing wrong it would be appreciated.

[Torsten]

You did nothing wrong, but why using the inexact gradient function if all derivatives are available from the ode calculation ?

[Paul Miller]

I am just unfamiliar using ODE and didn't know how to access the derivatives, and since it worked the first two times using gradient I was not sure why it didn't work for this equation.

[Torsten]

and since it worked the first two times using gradient I was not sure why it didn't work for this equation.

As said, it also works for this equation, but the results are not as exact as the direct derivatives from the ODE solution.

Answer 2

[John D'Errico]

Reemmber that the gradient function usues an APPROXIMATION. And, that at the ends of the series, that approximation is less good than it is in the middle part.

Worse, differentiation is a noise amplifying process. Tiny errors in the data will be amplified. Then taking the second derivative, and you get amplification squared.

Yes, gradient will sometimes work. Sometimes you get lucky. On a bad day, with a problem that is not quite so easy to work with, and where the approximations used by gradient are less accurate, you will see problems. That should be no surprise.

So what did you see? Crap happens with gradient near the endpoints. Again, gradient does NOT compute the exact derivative, just an approximation.

The point is, you got exactly what you might expect (ok, what you should have expected.) A far better estimate of the derivatives will be gained from the ODE function itself, in your case, f3. This is because that is the function ode15s used to solve the ODE itself.

[Paul Miller]

I see, so using @Torsten's method of plugging the values back into f3 minimizes the error which I introduced using the gradient function, and it only worked well when I did it before because I was using it for the right functions. Thank you.