LuisRei.com

I finally decided to take the Numerical Methods (MNUM) course. It turns out it’s a lot more fun than I thought. There is programming involved but you can chose to use whatever language you want. This is yet another nice excuse for me to use Python instead of C++ or Java. Last semester I was able to use Python to implement the game logic for Software Application Laboratory (LAS), which is mostly an OpenGL course with IPC via sockets thrown into the mix, and to write an article on dynamic languages (focusing mostly on Python) for Software Engineering (ESOF).

But back to this semester, 3 classes into the semester and the teacher is already said something like “I’m going to learn python now. I didn’t believe when I heard someone saying it was the best language in the world, but now I see there might be some truth to that claim”. That and I suspect his next laptop might be a macbook but that’s another story.

There are a few things that make Python great for Numerical Methods. In my opinion, Python’s clear, easy to understand, syntax is the most important one.It makes algorithms easier to implement. The syntax ends up being very close to language neutral pseudocode available in numerical methods books. Also Python’s datatypes as well as those provided by other libraries can be very useful.

The following code implements the stuff in chapter 2 (determining zeros) of the course. The methods implemented are Bisection, Rope and Newton. The function returns both the solution and the number of iterations necessary to get to that solution.

UPDATE: forgot the book - Numerical Methods in Engineering with Python

Appendix A - mnum2.py

from math import log

def bisect(f, a, b, e):
	""" Determines zero between a and b using Bisection. """
	n = 0
	fa = f(a)
	if fa == 0.0: return (a, n)
	fb = f(b)
	if fb == 0.0: return (b, n)

	while (abs(a-b) > e):
		c = 0.5*(a+b)
		fc = f(c)

		if fc == 0.0: return (c, n)
		n = n + 1
		if fb*fc < 0.0:
			a = c
			fa = fc

		else:
			b = c
			fb = fc

	if fa < fb:
		return (a, n)
	else:
		return (b, n)

def rope(f, a, b, e):
	""" Determines zero between a and b using the Rope methode. """
	n = 0
	fa = f(a)
	if fa == 0.0: return (a, n)
	fb = f(b)
	if fb == 0.0: return (b, n)

	while (abs(a-b) > e):
		c = (a*fb - b*fa) / (fb - fa)
		fc = f(c)
		if fc == 0.0: return (c, n)
		n = n + 1
		if fb*fc < 0:
			a = c
			fa = fc

		else:
			b = c
			fb = fc

	if fa < fb:
		return (a, n)
	else:
		return (b, n)

# Note: must verify that for the function f and guess c
#		the method will _converge_.
def newton(f, df, c, t):
	""" Determines zero between a and b using Newton """
	n = 0
	fc = f(c)
	if fc == 0.0: return (c, n)

	while (True):
		fc = f(c)
		dfc = df(c)
		if dfc == 0:
			print "dfc is 0"
			return (0, -1)

		dc = -fc/dfc

		c = c + dc
		n = n + 1
		if abs(dc) < t: return (c, n)

##Tests
#def f(x): return -log(x)+4.0
#def df(x): return -1.0/x
#x= bisect(f, 1, 70, 0.00000001)
#print x
#x = rope(f, 1, 70, 0.00000001)
#print x
#x = newton(f, df, 0.1, 0.0001)
#print x