#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. # #This is the basic version, using no plotting and using only #Python language constructs that are more or less common to #all programming languages. import math #Define the function you are integrating. To #keep things simple, this function has no parameters. #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 #Do the loop to find the solution starting from the #initial condition x = xstart y = ystart print "Euler Method" print "x y yExact y-yExact" for i in range(nsteps): slope = F(x,y) x = x + dx y = y + slope*dx yExact = ystart*math.exp(-x**2/2.) print x,y,yExact,y-yExact #Do a Midpoint method integration and save results for plotting #Do the loop to find the solution starting from the #initial condition. For dx = .5, it looks like the #midpoint method is actually less accurate than the Euler method, #out around x = 5. This is because the midpoint method has an #instability for this particular function. See the tutorial notes #for details. The instability problem goes away if you take a #smaller dx. x = xstart y = ystart print "Midpoint Method" print "x y yExact y-yExact" for i in range(nsteps): #First guess slope = F(x,y) ymid = y + .5*slope*dx #Now correct the slope slopeMid = F(x+dx/2.,ymid) #Now do the increment y = y + slopeMid*dx x = x+dx yExact = ystart*math.exp(-x**2/2.) print x,y,yExact,y-yExact #Do a Runge-Kutta method integration and save results for plotting #Do the loop to find the solution starting from the #initial condition x = xstart y = ystart print "Runge-Kutta Method" print "x y yExact y-yExact" for i in range(nsteps): #Define some fractional steps h = dx hh=h*0.5 h6=h/6.0 xh=x+hh #First slope estimate dydx = F(x,y) yt = y+hh*dydx #Next slope estimates dyt = F(xh,yt) yt =y+hh*dyt dym = F(xh,yt) yt =y+h*dym dym = dym + dyt dyt = F(x+h,yt) #Preliminaries done. Now combine slope estimates to get #the final increment y = y+ h6*(dydx+dyt+2.0*dym) x = x+h # yExact = ystart*math.exp(-x**2/2.) print x,y,yExact,y-yExact