0

I try to compare the analytical derivative of a sine function and approximative derivative using backward,forward,central,3.order backward difference method.I have huge difference between them. It must be mistake. Could you give me a hint where I made the mistake? Thank you in advance!

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt

from math import pi
from numpy import *



# spatial domain
xmin = 0
xmax = 1
n = 50 # num of grid points

# x grid of n points
x = np.linspace(xmin, xmax, n+1);

k=2


def f1(x):
    return np.sin(2*pi*k*x)
def f2(x):
    return np.pi*k*2*cos(2*pi*k*x)   #analytical derivative of f1 function

f = zeros(n+1,dtype=float)        # array to store values of f
fd=zeros(n+1,dtype=float)         # array to store values of fd
dfrwd = zeros(n+1,dtype=float)    # array to store values of calulated derivative with forward difference
dbwd = zeros(n+1,dtype=float)     # array to store values of calulated derivative with backward difference
dzd = zeros(n+1,dtype=float)      # array to store values of calulated derivative with central difference
ddo = zeros(n+1,dtype=float)      # array to store values of calulated derivative with 3.order backward difference

for i in range(0,n): # adds values to arrays for x and f(x)

    f[i] = f1(x[i])
    fd[i] = f2(x[i])
step=f[2]-f[1]

# periodic boundary conditions


dfrwd[n] = (f[0]-f[n])/step

dbwd[0] = (f[0]-f[n])/step

dzd[0]=(f[1]-f[n])/(2*step)
dzd[n]=(f[0]-f[n-1])/(2*step)

ddo[n] = (2*f[0]+3*f[n]-6*f[n-1]+f[n-2])/(6*step)
ddo[1] = (2*f[2]+3*f[1]-6*f[0]+f[n])/(6*step)
ddo[0] = (2*f[1]+3*f[0]-6*f[n]+f[n-1])/(6*step)


for i in range(0,n-1): # add values to array for derivative with forward difference
    dfrwd[i] = (f[i+1]-f[i])/step
for i in range(1,n): # add values to array for derivative with backward difference
    dbwd[i] = (f[i]-f[i-1])/step
for i in range(1,n-1): # add values to array for derivative with central difference
    dzd[i] = (f[i+1]-f[i-1])/(2*step)
for i in range(2,n-1): # add values to array for derivative with 3.order backward difference
    ddo[i] = (2*f[i+1]+3*f[i]-6*f[i-1]+f[i-2])/(6*step)




plt.plot(x,fd)
plt.plot(x,dzd)
alimuradpasa
  • 135
  • 1
  • 10

1 Answers1

1

You might want to compute step as x[1]-x[0].


Clean up your imports. You import and use numpy as np but also as a general import of all content, np.sin in f1 but only cos in f2. You use pi from math and np.pi in the same formula. This all is not wrong, but confusing and makes debugging harder.

You can use numpy vector operations, f = f1(x); fd = f2(x). Likewise the numerical derivatives can be computed via index shifts. Or you could implement them as dfrwd = (f1(x+step)-f1(x))/step etc. A little less efficient, but cleaner to read.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51