Numerical Analysis

[jamd200]

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)

[GDICommander]

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

[jamd200]

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.

[jamd200]

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

[GDICommander]

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.

[jamd200]

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