Forward, Central and Extrapolated Differences

mail@pastecode.io avatar
unknown
plain_text
2 years ago
2.2 kB
2
Indexable
Never
"""Numerical Differentiation"""

def FD(f, x, h): #forward method
    return (f(x + h) - f(x))/h

def CD(f, x, h): #forward method
    return (f(x + (0.5)*h) - f(x-(0.5)*h))/(h)

def ED(f, x, h): #forward method
    return (1/3)*(4*(CD(f,x,(0.5*h)))- CD(f,x,h))


from numpy import cos, sin, exp, logspace
from matplotlib.pyplot import loglog, xlabel, ylabel, title, legend
%matplotlib inline

test_f = np.cos   # choose an appropriate function here
x0 = 0.1       # where will we evaluate this function?
#h = 0.01 #set value of h
fx0p = (-3.19848111*(10**-6)) #for a value of h = 0.01, the true value of the derivative of test_f at x0

hh = logspace(-1, -17, 17) # same syntax as linspace: this gives us a range from 10^-1 to 10^-17 with 17 points.

fd_errors = []
cd_errors = []
ed_errors = []

# We will collect the epsilon values for the FD method in this list. 
# You may like to set up similar lists for other methods.

for h in hh:
    fd_exact = (-4.5692555*(10**-5))
    fd_estimate = FD(np.cos, 0.1, 10**(-1))
    cd_exact = (-3.04617226*(10**-5))
    cd_estimate = CD(np.cos, 0.1, 10**(-1))
    ed_exact = (-2.03078228*(10**-5))
    ed_estimate = ED(np.cos, 0.1, 10**(-1))

    fd_error = (fd_estimate - fd_exact)/(fd_exact)
    cd_error = (cd_estimate - cd_exact)/(cd_exact)
    ed_error = (ed_estimate - ed_exact)/(ed_exact)

    print("{:5.0e} {:10.9f} {:10.9f}".format(h, fd_estimate, fd_error))
    #print("{:5.0e} {:10.9f} {:10.9f}".format(h, cd_estimate, cd_error))
    #print("{:5.0e} {:10.9f} {:10.9f}".format(h, ed_estimate, ed_error))

    fd_errors.append(abs(fd_error)) #the absolute value takes the modulus so makes negative values positive values
    cd_errors.append(abs(cd_error))
    ed_errors.append(abs(ed_error))

loglog(hh, fd_errors, 'o-', label="FD") # Same syntax as plot. The label is used in the legend.
loglog(hh, cd_errors, '*-', label="CD")
loglog(hh, ed_errors, 'H-', label="ED")

xlabel('step size $h$') # Note that we can include LaTeX-style maths within dollar signs.
ylabel('absolute error $|\epsilon|$')

title('Observing the absolute error for varying step-size.') # Include an appropriate string for a graph title here.
legend() #key