#This scripts illustrates how to use #the integrator in ClimateUtilities to #numerically solve the ordinary differential equation # # dY/dx = f(x,Y) # #for the function Y(x) import math from ClimateUtilities import * #Define the function that computes the derivatives. #the third argument allows us to pass parameters to #the function. Here, the parameter is a real number, #but it can be any python object at all (e.g. a list of #numbers, or even a function). def f(x,Y,a): return -a*x*Y #Create an instance of the integrator object dx = .1 #The increment to use. #Smaller values give more accurate results Ystart = 1. #The initial value of Y xstart = 0. #The initial value of x # m = integrator(f,xstart,Ystart,dx) # #Set the parameter value a = 2. m.setParams(a) #Sets the parameter to the value of a # #Make lists to save the results for plotting xList = [xstart] YList = [Ystart] # #Step through until x exceeds the desired target x = xstart while x < 5.: x,Y = m.next() # next() can take a new value of dx as optional argument #Put the values in the results list xList.append(x) YList.append(Y) #Compute the exact solution for comparison YexactList = [Ystart*math.exp(-a*x*x/2.) for x in xList] #Put results in a Curve object and plot c = Curve() c.addCurve(xList,'x') c.addCurve(YList,'Y','Computed') c.addCurve(YexactList,'Yexact','Exact') w = plot(c) #With dx = .1, the approximate solution is so accurate you #can't see the difference between the exact and approximate #curves. Try larger values of dx to see how the accuracy #degrades at larger dx. The error shows up first at the #larger values of x #Alternately, we can compute the difference and plot that: c1 = Curve() c1.addCurve(xList,'x') c1['diff'] = c['Y']-c['Yexact'] #creates new variable named 'diff', #which is the difference between #Y and Yexact w1 = plot(c1)