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)

## All 5 Replies

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 meeting, learning, and sharing knowledge.