#Tutorial on numerical integration of first order ODE. # #We integrate the first order ODE # dy/dx = -x*y #numerically using the Euler method, the midpoint method, #and the Runge-Kutta method, and compare with the exact #solution. # # #**ToDo: Add a section at the end showing how to pass a parameter #as a global. #**ToDo: Write a stripped down version using just basic loops and no graphics. import math from ClimateUtilities import * #============================================================= #This class defines an object which implements numerical integration #of a first order ordinary differential equation #of the form dy/dx = F(x,y). x is the independent variable and #y is the dependent variable. # #See the main script below for an example of how the class is used. #This class is a stripped-down version of the Runge-Kutta integrator class #called "integrator," which is defined in ClimateUtilities. #The first argument of the function "derivs" is the independent #variable and the second argument is the dependent variable # class ODEsolve: def __init__(self, derivs,xstart,ystart,dx=None): self.derivs = derivs self.x = xstart self.y = ystart self.dx = dx #Method for integration using midpoint algoritm #Note that the increment dx is an optional second argument #Method for integration using Euler algorithm def nextEuler(self,dx = None): if not (dx == None): self.dx = dx h = self.dx dydx = self.derivs(self.x,self.y) self.y =self.y + h*dydx self.x = self.x + h return self.x,self.y #Method for integration using the midpoint algorithm def nextMidpoint(self,dx = None): if not (dx == None): self.dx = dx h = self.dx dydx = self.derivs(self.x,self.y) ymid =self.y + .5*h*dydx self.y = self.y + h*self.derivs(self.x+ .5*h,ymid) self.x = self.x + h return self.x,self.y #Method for integration using Runge-Kutta def nextRK(self,dx=None): if not (dx == None): self.dx = dx h = self.dx hh=h*0.5; h6=h/6.0; xh=self.x+hh; dydx = self.derivs(self.x,self.y) yt = self.y+hh*dydx dyt = self.derivs(xh,yt) yt =self.y+hh*dyt dym = self.derivs(xh,yt) yt =self.y+h*dym dym += dyt dyt = self.derivs(self.x+h,yt) self.y += h6*(dydx+dyt+2.0*dym) self.x += h return self.x,self.y #============================================================= #-----Main Script starts here---- #Define the function you are integrating. To #keep things simple, this function has no parameters. #Later, we'll learn how to pass extra parameters #to the function. The analytic solution to this #problem is y(x) = y(0)*exp(-x**2/2) def F(x,y): return -x*y #Initial conditions xstart = 0 ystart = 1. # #Increment. Try re-running the script with different #values dx = .5 # #How many steps? nsteps = 10 #Try 15 steps, to show instability of midpoint method #Do an Euler integration and save results for plotting xList = [xstart] #Same xList for all cases yEulerList = [ystart] mEuler = ODEsolve(F,xstart,ystart,dx) #Create an instance of an integrator #Do the loop to find the solution starting from the #initial condition for i in range(nsteps): x,y = mEuler.nextEuler() xList.append(x) yEulerList.append(y) #Now re-do it for the midpoint method yMidpointList = [ystart] mMidpoint = ODEsolve(F,xstart,ystart,dx) #Create an instance of an integrator #Do the loop to find the solution starting from the #initial condition for i in range(nsteps): x,y = mMidpoint.nextMidpoint() yMidpointList.append(y) #Now do it using the more accurate Runge-Kutta integrator yRKList = [ystart] mRK = ODEsolve(F,xstart,ystart,dx) #Do the loop for i in range(nsteps): x,y = mRK.nextRK() yRKList.append(y) #Now make a list of the analytic solution yExact = [math.exp(-x*x/2.) for x in xList] #Put it in a Curve object for plotting and output c = Curve() c.addCurve(xList,'x') c.addCurve(yEulerList,'Euler') c.addCurve(yMidpointList,'Midpoint') c.addCurve(yRKList,'RK') c.addCurve(yExact,'Exact') # #Annotations c.Xlabel = 'x' c.Ylabel = 'y' c.PlotTitle = 'Approximate and exact solutions for dy/dx = -xy' # #Plot it, and then write it out for use in an external #plotting program plot(c) #Now, it's a bit hard to see the difference between the curves, #so instead, let's plot the difference between the approximate #and exact solutions. This calculation uses the fact that the #Curve object stores the data in a way that can be used to #do array arithmetic. c1 = Curve() c1.addCurve(xList,'x') c1.addCurve(c['Euler']-c['Exact'],'EulerDiff') c1.addCurve(c['Midpoint']-c['Exact'],'MidpointDiff') c1.addCurve(c['RK']-c['Exact'],'RKDiff') #Annotations c1.Xlabel = 'x' c1.Ylabel = 'y - yExact' c1.PlotTitle = 'Approximate minus exact solutions for dy/dx = -xy' plot(c1)