3.2. Defining your own harmonic oscillator equations
Up to now, we have loaded the predefined harmonic oscillator equations with the line
from pyoomph.equations.harmonic_oscillator import HarmonicOscillator
However, one of the main features of pyoomph is actually the definition of arbitrary equations. Instead of including the predefined equation class, we will write our own class, inheriting from the generic ODEEquations. This will serve as an example of how to express arbitrary ODEs within pyoomph. At the same time, let us generalize the equation to include damping \(\delta\) and driving \(f(t)\) as follows
The corresponding code reads as follows:
from pyoomph import * # Import pyoomph
from pyoomph.expressions import * # Import some additional things to express e.g. partial_t
# We define a new class called HarmonicOscillator, which is inherited from the generic ODEEquations
class HarmonicOscillator(ODEEquations):
# Constructor, allows us to set some parameters like the name of the variable, omega, damping and driving
def __init__(self,*,name="y",omega=1,damping=0,driving=0):
super(HarmonicOscillator,self).__init__()
self.name=name #Store these as members of the equation object
self.omega=omega
self.damping=damping
self.driving=driving
# This function is called to define all fields in this ODE (system)
def define_fields(self):
self.define_ode_variable(self.name)
# This function will finally define the equations
def define_residuals(self):
y=var(self.name) # bind the local variable y to the ODE variable
# Write the equation in residual form, i.e. lhs-rhs=0
residual=partial_t(y,2)+2*self.damping*partial_t(y)+self.omega**2 *y - self.driving
# And add the residual to the equation. Here, we have to project it on the test function.
self.add_residual(residual*testfunction(y))
In the __init__() method, we take optional keyword arguments which have default values. These are then stored as members of the class. Of course, we have to call again the constructor of the super-class ODEEquations using the Python builtin super.
In the next step, the method define_fields() is overloaded. Inside this function, all unknowns of the ODE system have to be defined. Here, it is just the single unknown \(y\).
Finally, the nitty-gritty, the equation has to be defined. This happens in the method define_residuals(). To that end, the equation first has to be cast to a residual form, which can be done by putting all terms on one side of the equation:
Of course, only the lhs of this equation is of relevance. Before adding it to the system of equations with the add_residual() method, it must be multiplied by a so-called test function, which can be obtained by the function testfunction(). This all happens in the line self.add_residual(residual*testfunction(y)). Here, we only have one equation to be solved, i.e. the harmonic oscillator. In general, one will have a system of ODEs, i.e. multiple equations for multiple unknowns. This is where test functions become important, which will be explained in the next section.
The remainder of the code is very similar to before, but now also damping and driving will be considered:
# The remainder is almost the same is in the example nondim_harmonic_osci.py
class HarmonicOscillatorProblem(Problem):
def __init__(self):
super(HarmonicOscillatorProblem,self).__init__()
self.omega=1
self.damping=0.1 #But we add some default damping here
t=var("time")
self.driving=0.2*cos(0.2*t) #and some driving
def define_problem(self):
eqs=HarmonicOscillator(omega=self.omega,damping=self.damping,driving=self.driving,name="y") #We also pass the damping and driving here
eqs+=InitialCondition(y=1-var("time"))
eqs+=ODEFileOutput()
self.add_equations(eqs@"harmonic_oscillator")
if __name__=="__main__":
with HarmonicOscillatorProblem() as problem:
problem.run(endtime=100,numouts=1000)
The output is plotted in Fig. 3.2.
Fig. 3.2 Output for the user-defined harmonic oscillator equation with damping and driving.