Dfenestrator

By the numbers

That's the last straw. Sage can't handle fractional exponents on negative numbers? Pfaw! No big, can just do it in Mathematica. Functional is more my style, anyway. Dynamically-typed imperative languages are for the lazy or spastic.

So I bequeath to you the Python code of my aborted Sage notebook for my numerical differential equation solver. Euler's method only. Was going to do Improved Euler's and RKx but I'll save it for Mathematica.

Use it together.

Use it in peace.

Cell 1

def Euler(f, x0, y0, xf, h):
    
    x  = float(x0)
    y  = float(y0)
    xf = float(xf)
    h  = float(h)
    
    n = int((xf - x)/h)
        
    xs = [x]
    ys = [y]
    
    for i in range(0, n):
        y = y + h * f(x,y).subs(x=x, y=y).n()
        x = x + h
        xs.append(x)
        ys.append(y)
        
    return [xs,ys]

Cell 2

def TabList(y, headers = None):
    '''
    Converts a list into an html table with borders.
    '''
    s = '<table border = 1>'
    if headers:
        for q in headers:
            s = s + '<th>' + str(q) + '</th>'
    for x in y:
        s = s + '<tr>'
        for q in x:
            s = s + '<td>' + str(q) + '</td>'
        s = s + '</tr>'
    s = s + '</table>'
    return s

Cell 3

class NumericalSolution:
    def __init__(self, fp, x0, y0, xf, h, color = (0,0,1)):
        self.fp = fp
        self.x0 = x0
        self.y0 = y0
        self.xf = xf
        self.h = h
        self.color = color
    def Solve(self):
        raise MethodNotImplemented()
    def ToLine(self):
        X = 0; Y = 1;
        result = self.Solve()
        return line( [[result[X][i], result[Y][i]] for i in range(len(result[X]))], rgbcolor=self.color )

class MethodNotImplemented: pass

Cell 4

class EulerNumericalSolution(NumericalSolution):
    def __init__(self, fp, x0, y0, xf, h, color = (0,0,1)):
        NumericalSolution.__init__(self, fp, x0, y0, xf, h, color)
    def Solve(self):
        #return Euler( self.fp, self.x0, self.y0, self.xf, self.h)
        return Euler( (lambda x,y: self.fp), self.x0, self.y0, self.xf, self.h)

Cell 5

SolverTypes = ["Euler"] #, "Improved Euler", "RK2"]

def GenerateNumericalSolution( fp, solver, color = (0,0,1) ):        
    if solver == "Euler": 
        return EulerNumericalSolution(fp, x0, y0, xf, h, color)
    else:
        raise InvalidSolver()
        
class InvalidSolver: pass

Cell 6

def GetNumericalSolutionLines():
    lines = solutions[0].ToLine()
    for i in range(1, len(solutions)):
        lines += solutions[i].ToLine()
    return lines

Cell 7

f = 0
fColor = (1,0,0)
@interact
def DefineExactFunction( f_in = input_box(-cos(x)+1.0, label = 'f(x) = '),
                         color = color_selector((1,0,0), label='', widget='farbtastic', hide_box=False) ):        
    global f
    global fColor
    f = f_in
    fColor = color
    return

Cell 8

x0 = 0.0
y0 = 0.0
xf = 10.0
h  = 0.1

@interact
def _( _x0 = input_box( x0, label = 'x<sub>0</sub> '),
       _y0 = input_box( y0, label = 'y<sub>0</sub> '),
       _xf = input_box( xf, label = 'x<sub>f</sub> '),
       _h  = input_box(  h, label = 'h ') ):
    global x0; global y0; global xf; global h    
    x0  = _x0; y0  = _y0; xf  = _xf; h = _h;
    
    return

Cell 9

#x0 = 0.0
#y0 = 0.0
#xf = 10.0
#h  = 0.1
#
#@interact
#def _( _x0 = input_box( x0, label = 'x<sub>0</sub> '),
#       _y0 = input_box( y0, label = 'y<sub>0</sub> '),
#       _xf = input_box( xf, label = 'x<sub>f</sub> '),
#       _n  = input_box( int(xf-x0/h), label = 'number of steps ') ):
#    global x0; global y0; global xf; global h    
#    x0  = _x0; y0  = _y0; xf  = _xf; h = (_xf - _x0) / _n;
#    
#    return

Cell 10

# *** Enter Numerical Solutions ***
solutions = []
var('x,y')
@interact
def _(
       fp     = input_box(sin(x), label = "f'(x,y)="),
       solver = selector(SolverTypes, label="Solver="),
       color  = color_selector((0,0,1), label='', widget='farbtastic', hide_box=False),
       action = selector(["Edit","Add","Clear All"], label='', buttons=True)
     ):
    global solutions
    if action == "Add":
        solutions.append( GenerateNumericalSolution( fp, solver, color ) )
    elif action == "Clear All":
        solutions = []
    else:
        pass

    for i in range(len(solutions)):
        html('<font color=\"%s\">$f\'\;(x)\;=\;%s,\;h=%f$</font>' % (solutions[i].color.html_color(), latex(solutions[i].fp), solutions[i].h))
        
    return

Cell 11

@interact
def _( showExact = checkbox (True, "Show exact") ):
    X = 0
    Y = 1
    #try:
    if showExact:
        show( plot(f,x0,xf,rgbcolor=fColor) + GetNumericalSolutionLines() )
    else:
        show( GetNumericalSolutionLines() ) 
    if showExact:
        html('<font color=\"%s\"><center>$f\'\;(x)\;=\;%s\;%s$</center></font>' % ( fColor.html_color(), latex(f), "(exact)" ))    
    for i in range(len(solutions)):
        solution = solutions[i].Solve()
        tableRange = range(len(solution[X]))
        html('<font color=\"%s\"><center>$f\'\;(x)\;=\;%s,\;h=%f$</center></font>' % (solutions[i].color.html_color(), latex(solutions[i].fp), solutions[i].h))    
    #except:
    #print "< < No Functions Configured > >"

Cell 12

@interact
def _( showExact = checkbox (True, "Show exact") ):
    if len(solutions) == 0:
        print "< < No Functions Configured > >"
    else:
        X = 0
        Y = 1
        for i in range(len(solutions)):
            solution = solutions[i].Solve()
            tableRange = range(len(solution[X]))
            html('<font color=\"%s\">$f\'\;(x)\;=\;%s,\;h=%f$</font>' % (solutions[i].color.html_color(), latex(solutions[i].fp), solutions[i].h))
            html('')
            if showExact:
                html(TabList([[j, solution[X][j], solution[Y][j], f(x=solution[X][j]).n()] for j in tableRange], headers = ['step','x','y', 'y (exact)']))
            else:
                html(TabList([[j, solution[X][j], solution[Y][j]] for j in tableRange], headers = ['step','x','y']))
            html('')

Cell 13

@interact
def _( showExact = checkbox (True, "Show exact"),
       displayInterval = input_box( 0.1, label = 'Interval ') ):
    if len(solutions) == 0:
        print "< < No Functions Configured > >"
    else:
        X = 0
        Y = 1

        headers = ['x']
        for i in range(len(solutions)):
            yHeader = "y (h=%f)" % solutions[i].h
            headers.append(yHeader)
        if showExact:
            headers.append('y (exact)')
        
        solution = []
        decimatedSolution = []
        for i in range(len(solutions)):
            solution.append(solutions[i].Solve())
            decimatedSolution.append([])
        #if showExact:
            #solution.append([])
        
        numIntervals = int((xf - x0)/displayInterval)
        
        for i in range(len(solution)):  
            numPointsInInterval = len(solution[i][Y])/numIntervals
            for j in range(numIntervals+1):
                decimatedSolution[i].append(solution[i][Y][int(j*numPointsInInterval)])
        
        vals = []
        for i in range(len(decimatedSolution[0])): # range(len(solution[0][X])-1):
            slice = []
            slice.append(i*displayInterval)
            for j in range(len(solution)):
                slice.append(decimatedSolution[j][i])
            if showExact:
                slice.append(f(x=i*displayInterval).n())
            vals.append(slice)
        html(TabList(vals, headers))

-- Posted 12:03:00 AM by Paul Hrissikopoulos