from ClimateUtilities import *
import math
#Set up a function that computes the slope
#as a function of the dependent and independent
#variables.
#
#Here, the proportionality constant "a" is a parameter.
#This version shows how to pass a parameter using
#as an argument to the slope function. This is convenient
#if you ever want to put your solver inside another function
#and pass the parameter as an argument to that function.
#(Doing that via globals is awkward). Note that the
#parameter can actually be any object. Here it is just
#a real number but below we'll show an example where
#an object is used to pass multiple values (even functions)
#to the slope function.
#Here Y is the dependent variable and t is the independent variable
#The independent variable comes first in the argument list
def dYdt(t,Y,a):
return -a*Y
#Now create an instance of the integrator for the function dYdt
tStart = 0. #Starting value of independent variable
YStart = 1. #Starting value of dependent variable
dt = .1 #The increment to use. Smaller is more accurate
m = integrator(dYdt,tStart,YStart,dt)
#Set the parameter value. (You could skip this and instead
#write the function dY/dt without the parameter in the argument
#list, and treat the parameter as a global instead)
a = 2.
m.setParams(a)
#The following loop carries out the integration to t = 10.,
#and just prints out the results together with the analytic
#solution for comparison. I use a "for" loop here for simplicity,
#but a better way would be to use a "while" loop to test when
#to stop.
nsteps = int(10./dt)
for i in range(nsteps):
t,Y = m.next() #The next() method advances the values by one time step
print t,Y,math.exp(-a*t)
#----------------Example 2: Using an object to pass multiple parameters.
#
#Sometimes, the computation of the slope function depends on many
#different constants, and you might even need to specify a function
#(e.g. the time dependence of something). You could do this by
#making the parameter a list (which could be a list of any kind of
#objects) and refer to the various elements as a[0], a[1], etc. in
#your slope function. However, here's a nicer way that allows you
#to refer to the various things by name. Easier to keep track of
#what's what that way.
#First we set up a dummy object, with no members. The class Dummy()
#is defined in ClimateUtilities
params = Dummy()
#Now let's say we want to pass a function f(t) to our slope function
#First define it
def f(t):
return math.sin(t)
#Now create members of the object "params". You do this by
#just assigning values! This can be done at any point before
#you make use of the integrator object, and parameters can
#be changed at any time in the course of the integration
params.TimeFunction = f
#Suppose we also want to pass a real number "a"
params.a = 2.
#Now, this is how you write the slope function to use the params object
#Note that it is a dummy variable, so it can be called anything. It
#does not have to have the same name as the object we just defined!
def dYdt(t,Y,InParams):
a = InParams.a
f = InParams.TimeFunction
return -a*f(t)*Y
#Now set up the integration as before
#Now create an instance of the integrator for the function dYdt
tStart = 0. #Starting value of independent variable
YStart = 1. #Starting value of dependent variable
dt = .1 #The increment to use. Smaller is more accurate
m = integrator(dYdt,tStart,YStart,dt)
#The parameter value was already set above, but we need
#to tell the integrator to use it
m.setParams(params)
# Note that if you reset some member of params, you don't
#need to do a setParams again. e.g. if you want to
#change the value of a to 3., just do params.a = 3.
#The following loop carries out the integration to t = 10.,
#and just prints out the results together with the analytic
#solution for comparison. I use a "for" loop here for simplicity,
#but a better way would be to use a "while" loop to test when
#to stop.
nsteps = int(10./dt)
for i in range(nsteps):
t,Y = m.next() #The next() method advances the values by one time step
print t,Y