So far, my simple numerical analysis program has differentiation, integration, first-order ODEs and Taylor Series. What else can I add to this?

from math import*

"""Single-variable calculus."""

#DIFFERENTIATION

def gen_differentiation(f,x,h):
    #finds the first derivative of any mathematical function
    return (f(x+h)-f(x))/(h)

def special_differentiation(f,x):
    #finds the first derivative of specific functions
    if (f == sin(x)):
        return cos(x)
    elif (f == cos(x)):
        return -sin(x)
    elif (f == tan(x)):
        return 1/(sin(x)*sin(x))
    elif (f == exp(x)):
        return exp(x)
    elif (f == log(x,e)):
        return 1/x
    else:
        return "Not a special function."

def second_derivative(f,x,h):
    #Finds the second derivative of any mathematical function
    return (f(x+h)-2*f(x)+f(x-h))/(h*h)

    
def third_derivative(f,x,h):
    #Finds the third derivative of any mathematical function
    numer = (h*h)*(gen_differentiation(f,x+h,h)-2*gen_differentiation(f,x,h)+\
            gen_differentiation(f,x-h,h))
    denom = h**4
    return numer/denom

#INTEGRATION
    
def integration(f,a,b): #Using Simpson's Rule
    return ((b-a)/6)*(f(a)+4*f(0.5*(a+b))+f(b))

def integration1(f,a,b): #Using Trapezoidal Rule (general)
    return (b-a)*0.5*(f(a)+f(b))

def integration2(f,a,b,n): #Using Trapezoidal Rule (strips)
    clist = [f(a+i*((b-a)/n)) for i in range(1,n)]
    return ((b-a)/n)*(0.5*(f(a)+f(b))+ sum(clist))


#FIRST-ORDER DIFFERENTIAL EQUATIONS


def direct_method(dydx,x1,x2):
    #Only works with single-variable
    y = integration(dydx, x1, x2)
    return y


def midpoint(f,x0,y0,h,limit): #The midpoint method
    """f is the value of dy/dx."""
    while x0 <= limit:
        k1 = h*f(x0,y0)
        k2 = h*f(x0 + 0.5*h, y0 + 0.5*k1)
        y0 = y0 + k2
        x0 = x0 + h

    return y0

def runge_kutta(f,x0,y0,h,limit): #The Runge-Kutta 4th-order Method
    """f is the value of dy/dx."""
    while x0 <= limit:
        k1 = h*f(x0,y0)
        k2 = h*f(x0 + 0.5*h, y0 + 0.5*k1)
        k3 = h*f(x0 + 0.5*h, y0 + 0.5*k2)
        k4 = h*f(x0 + h, y0 + k3)

        y0 = y0+(k1/6)+(k2/3)+(k3/3)+(k4/6)
        x0 = x0 + h

    return y0


#TAYLOR SERIES

def taylor_series(f,x,a): #Taylor series up to the fourth term
    d = gen_differentiation(f,a,0.0000001)
    sd = second_derivative(f,a,0.0000001)
    td = third_derivative(f,a,0.0001)
    second_term = d*(x-a)
    third_term = (sd*(x-a)*(x-a))/factorial(2)
    fourth_term = (td*(x-a)**3)/factorial(3)
    series = [f(a), second_term, third_term, fourth_term]
    return sum(series)

Recommended Answers

I would add what is necessary to add... Who will be using your code?

Jump to Post

Have you tested all your functions with unit tests? (You can easily write unit tests in Python, see http://docs.python.org/library/unittest.html)

If you haven't tested what you wrote, then you may find a lot …

Jump to Post

All 5 Replies

I would add what is necessary to add... Who will be using your code?

I would add what is necessary to add... Who will be using your code?

This is something I'm doing on the side to practice my programming and to prepare myself a little for upcoming Math topics (I'm doing both Computer Science and Math). Also, I'm into creating simple scientific modules myself, instead of downloading SciPy or NumPy.

I changed the Runge-Kutta and midpoint functions to return a list of x-y coordinate points
for x-values from x0 to limit.

def midpoint(f,x0,y0,h,limit): #The midpoint method
    """f is the value of dy/dx."""
    point_lst = []
    point_lst.append((x0,y0))
    while x0 < limit:
        k1 = h*f(x0,y0)
        k2 = h*f(x0 + 0.5*h, y0 + 0.5*k1)
        y0 = y0 + k2
        x0 = x0 + h
        point_lst.append((x0,y0))

    #Returns a list of x-y points for a graph 
    return point_lst

def runge_kutta(f,x0,y0,h,limit): #The Runge-Kutta 4th-order Method
    """f is the value of dy/dx."""
    point_lst = []
    point_lst.append((x0,y0))
    while x0 < limit:
        k1 = h*f(x0,y0)
        k2 = h*f(x0 + 0.5*h, y0 + 0.5*k1)
        k3 = h*f(x0 + 0.5*h, y0 + 0.5*k2)
        k4 = h*f(x0 + h, y0 + k3)

        y0 = y0+(k1/6)+(k2/3)+(k3/3)+(k4/6)
        x0 = x0 + h
        point_lst.append((x0,y0))

    #Returns a list of x-y points for a graph   
    return point_lst

Have you tested all your functions with unit tests? (You can easily write unit tests in Python, see http://docs.python.org/library/unittest.html)

If you haven't tested what you wrote, then you may find a lot of bugs inside your code when you'll be using it for real.

Thanks for the advice. I'm still a student, so I'm still learning the ropes.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts learning and sharing knowledge.