FOURTH DERIVATIVE - My graph is wrong, please let me know my mistake

Viewed 57

For the fourth derivative, this is my code:

import numpy as np
import matplotlib.pyplot as plt

pi=np.pi
x=np.arange(0,5*pi,0.03*pi)
y=np.cos(x)

h=x[4]-x[3]

"fourth derivative"

#first formula

der41=[]
x41=[]

for i in range(2,len(y)-2,1):
    der41.append((y[i+2]-4*y[i+1]+6*y[i]-4*y[i-1]+y[i-2])/(h**4))
    x41.append(x[i])

der41=np.array(der41)
x41=np.array(x41)

#Second formula

der42=[]
x42=[]

for i in range(3,len(y)-3,1):
    der42.append((-y[i+3]+12*y[i+2]+39*y[i+1]+56*y[i]-39*y[i-1]+12*y[i-2]+y[i-3])/(6*h**4))
    x42.append(x[i])

der42=np.array(der42)
x42=np.array(x42)

plt.plot(x,y,'g',x41,der41,'r',x42,der42,'b')
plt.title("Cuarta derivada")
plt.xlabel("Eje X")
plt.ylabel("Eje Y")
plt.grid()
plt.show()

This is the graph:

enter image description here

This code is based on the following formulas:

enter image description here

All three graphs should be the same "shape".

PD: Sorry if i have mistaken a word, please let me know. Thank you.

1 Answers

The printed formula that you show is incorrect. (I hope it's not out of a book.)

The coefficients should be (from wikipedia):

-y[i+3]+12*y[i+2]-39*y[i+1]+56*y[i]-39*y[i-1]+12*y[i-2]-y[i-3]

which gives (after adding a 0.1 offset to the curves to separate them):

enter image description here

Explanation:
More interesting, I think, is that for finite difference equations, all of the coefficients need to add to zero. Knowing this, you can see that the printed equation is incorrect at a glance. The other, and initial, key to this was noticing that you were plotting a cosine but that the magnitude was huge. What would make it so huge? Probably dividing a regular sized number (~1) by h**4, where h is small. So considering this, and, in fact, that this is always the situation in finite difference calculations, is where I realized that the sum of the coefficients must add to zero (which I state as a rule, but really I just made up, although I think it's probably a rule somewhere), which clearly isn't the case for the printed formula. Also, the 1D finite difference formulas are usually symmetric, and probably must be with the signs on each term if the coefficients are symmetric, so, that's another indication that the printed formula is incorrect. (I say all this assuming that you're learning these topics, and thought it would be useful to show it's just common sense type of reasoning, and nothing crazy or esoteric.)

Numpy:
Finally, I'll note that this will run faster, and be easier to write, take less memory, etc, if you keep it all in numpy, for example:

der42 = (-y[6:] + 12*y[5:-1] - 39*y[4:-2] + 56*y[3:-3] - 39*y[2:-4] + 12*y[1:-5] - y[:-6])/(6*h**4)
x42 = x[3:-3]
Related