1

Please look at my code and help me. I'm writing a code the approximates the Derivative of a function. To visually check how close my approximation is with the real derivative, I'm plotting both of them together.

My problem is when the function is not defined at zero, like 1/x with a derivative of '1/x^2.

Thanks in advance.

# -*- coding: utf-8 -*-
from pylab import *
import math
def Derivative(f,x, Tol = 10e-5, Max = 20):
try:
    k = 2
    D_old = (f(x+2.0**(-k))-f(x-2.0**-k))/(2.0**(1-k))
    k = 3
    D_new = (f(x+2.0**(-k))-f(x-2.0**-k))/(2.0**(1-k))

    E_old = abs(D_new - D_old)

    while True:
        D_old = D_new 

        k+=1

        D_new = (f(x+2.0**(-k))-f(x-2.0**-k))/(2.0**(1-k))
        E_new = abs(D_old - D_new)

        if E_new < Tol or E_new >= E_old or k >= Max:
            return D_old

except:
    return nan
def Fa(x):
     return math.sin(2*math.pi*x)
def Fap(x):
     return 2*math.pi*math.cos(2*math.pi*x)
def Fb(x):
    return x**2
def Fbp(x):
    return 2*x
def Fc(x):
    return 1.0/x
def Fcp(x):
    if abs(x)<0.01:
    return 0
    else:
        return -1.0/x**2
def Fd(x):
    return abs(x)
def Fdp(x):
    return 1 #since it's x/sqrt(x**2)

# Plot of Derivative Fa
xx = arange(-1, 1, 0.01)              # A Numpy vector of x-values
yy = [Derivative(Fa, x) for x in xx]  # Vector of f’ approximations
plot(xx, yy, 'r--', linewidth = 5)      # solid red line of width 5

yy2 = [Fap(x) for x in xx]
plot(xx, yy2, 'b--', linewidth = 2)     # solid blue line of width 2

# Plot of Derivative Fb

yy = [Derivative(Fb, x) for x in xx]  # Vector of f’ approximations
plot(xx, yy, 'g^', linewidth = 5)      # solid green line of width 5

yy2 = [Fbp(x) for x in xx]
plot(xx, yy2, 'y^', linewidth = 2)     # solid yellow line of width 2
Jomisilfe
  • 55
  • 1
  • 9
  • 1
    What is your question? – BrenBarn Mar 31 '14 at 17:49
  • My question is How to make a plot when the function is not defined at zero, like 1/x with a derivative of -1/x^2. I have tried to break the x-axis and it worked, but I want to know if there is another way to do it to make the code shorter! – Jomisilfe Mar 31 '14 at 17:58
  • Breaking the axis only makes the code one line longer: ```xx = xx[np.absolute(xx) > .001]``` – wwii Mar 31 '14 at 18:39

1 Answers1

1

"1/x" is an infinite function and you can not plot a function that is not defined to zero. You can only plot the function with a broken axis. For the broken axis, you can follow the suggestion of wwii in comments or you can follow this tutorial for matplotlib.

This tutorial show how you can just use two subplots to create the effect you desire.

Here the example code:

"""
Broken axis example, where the y-axis will have a portion cut out.
"""
import matplotlib.pylab as plt
import numpy as np


# 30 points between 0 0.2] originally made using np.random.rand(30)*.2
pts = np.array([ 0.015,  0.166,  0.133,  0.159,  0.041,  0.024,  0.195,
    0.039, 0.161,  0.018,  0.143,  0.056,  0.125,  0.096,  0.094, 0.051,
    0.043,  0.021,  0.138,  0.075,  0.109,  0.195,  0.05 , 0.074, 0.079,
    0.155,  0.02 ,  0.01 ,  0.061,  0.008])

# Now let's make two outlier points which are far away from everything. 
pts[[3,14]] += .8

# If we were to simply plot pts, we'd lose most of the interesting
# details due to the outliers. So let's 'break' or 'cut-out' the y-axis
# into two portions - use the top (ax) for the outliers, and the bottom
# (ax2) for the details of the majority of our data
f,(ax,ax2) = plt.subplots(2,1,sharex=True)

# plot the same data on both axes
ax.plot(pts)
ax2.plot(pts)

# zoom-in / limit the view to different portions of the data
ax.set_ylim(.78,1.) # outliers only
ax2.set_ylim(0,.22) # most of the data

# hide the spines between ax and ax2
ax.spines['bottom'].set_visible(False)
ax2.spines['top'].set_visible(False)
ax.xaxis.tick_top()
ax.tick_params(labeltop='off') # don't put tick labels at the top
ax2.xaxis.tick_bottom()

# This looks pretty good, and was fairly painless, but you can get that
# cut-out diagonal lines look with just a bit more work. The important
# thing to know here is that in axes coordinates, which are always
# between 0-1, spine endpoints are at these locations (0,0), (0,1),
# (1,0), and (1,1).  Thus, we just need to put the diagonals in the
# appropriate corners of each of our axes, and so long as we use the
# right transform and disable clipping.

d = .015 # how big to make the diagonal lines in axes coordinates
# arguments to pass plot, just so we don't keep repeating them
kwargs = dict(transform=ax.transAxes, color='k', clip_on=False)
ax.plot((-d,+d),(-d,+d), **kwargs)      # top-left diagonal
ax.plot((1-d,1+d),(-d,+d), **kwargs)    # top-right diagonal

kwargs.update(transform=ax2.transAxes)  # switch to the bottom axes
ax2.plot((-d,+d),(1-d,1+d), **kwargs)   # bottom-left diagonal
ax2.plot((1-d,1+d),(1-d,1+d), **kwargs) # bottom-right diagonal

# What's cool about this is that now if we vary the distance between
# ax and ax2 via f.subplots_adjust(hspace=...) or plt.subplot_tool(),
# the diagonal lines will move accordingly, and stay right at the tips
# of the spines they are 'breaking'

plt.show()
elviuz
  • 639
  • 1
  • 7
  • 26