3.4. Coupled harmonic oscillators
(3.4)\[\begin{split}\begin{split} \partial_t^2 y_1 + K_{11} y_1 + K_{12} y_2 &=0 \\ \partial_t^2 y_2 + K_{21} y_1 + K_{22} y_2 &=0 \end{split}\end{split}\]
Method I
We start by the naive way of implementing this, i.e. by defining a specific equation class, inherited from the generic ODEEquations class:
from pyoomph import * # Import pyoomph
from pyoomph.expressions import * # Import some additional things to express e.g. partial_t
class TwoCoupledHarmonicOscillator(ODEEquations):
def __init__(self,Kmatrix): #Required to pass the coupling matrix here
super(TwoCoupledHarmonicOscillator,self).__init__()
self.Kmatrix=Kmatrix #Store the matrix in the equations object
def define_fields(self):
self.define_ode_variable("y1") #define both unknowns
self.define_ode_variable("y2")
def define_residuals(self):
y1,y2=var(["y1","y2"]) #Bind both unknowns to variables
# Calculate the residuals
residual1=partial_t(y1,2)+self.Kmatrix[0][0]*y1+self.Kmatrix[0][1]*y2
residual2=partial_t(y2,2)+self.Kmatrix[1][0]*y1+self.Kmatrix[1][1]*y2
# Add both residuals
self.add_residual(residual1*testfunction(y1)+residual2*testfunction(y2))
# The remaining part is analogous to the normal oscillator
class TwoCoupledHarmonicOscillatorProblem(Problem):
def __init__(self):
super(TwoCoupledHarmonicOscillatorProblem,self).__init__()
self.Kmatrix=[[1,-0.5],[-0.2,0.4]] # Some coefficient matrix
def define_problem(self):
eqs=TwoCoupledHarmonicOscillator(self.Kmatrix)
eqs+=InitialCondition(y1=1,y2=0) #Setting initial condition
eqs+=ODEFileOutput()
self.add_equations(eqs@"coupled_oscillator")
if __name__=="__main__":
with TwoCoupledHarmonicOscillatorProblem() as problem:
problem.run(endtime=100,numouts=1000)
Here, in particular the line
self.add_residual(residual1*testfunction(y1)+residual2*testfunction(y2))
residual1 and residual2 in the code. The test functions \(V^{(j)}_i\) are just obtained by the calls of testfunction(). We actually solve now the residual of the first equation in (3.4) by adjusting \(y_1\) and the second equation by adjusting \(y_2\) until both equations are fulfilled.Method II
The previous method is straightforward, but would require to rewrite a lot of code if e.g. a third oscillator should be coupled to the system. Of course, one could use len to obtain the dimension of the coupling matrix Kmatrix and perform the calls for define_ode_variable() and the residual calculation in a loop. However, pyoomph also offers another way of coupling multiple equations. We have already seen that we can augment an equation e.g. with an InitialCondition object to set an initial condition or an ODEFileOutput object to make sure that output is written. By the very same way, also multiple equations can be coupled. To that end, we write an equation class for a general second order equation of the form
where \(T\) is a place holder for an arbitrary mathematical expression. The corresponding equation class looks like this:
from pyoomph import * # Import pyoomph
from pyoomph.expressions import * # Import some additional things to express e.g. partial_t
class SingleHarmonicOscillator(ODEEquations):
def __init__(self,name,terms): #Pass the name of the unknown and the terms T
super(SingleHarmonicOscillator,self).__init__()
self.name=name #Store the name of the unknown
self.terms=terms #and the terms to consider
def define_fields(self):
self.define_ode_variable(self.name)
def define_residuals(self):
y=var(self.name) #Bind the single variable
# Calculate the residuals
residual=partial_t(y,2)+self.terms #Just add the passed terms here
self.add_residual(residual*testfunction(y))
In the definition of the problem, we can now combine two instances of the SingleHarmonicOscillator class and by passing the correct names and additional terms, we can recreate the system (3.4). This is achieved by changing the method define_problem() of the TwoCoupledHarmonicOscillatorProblem as follows:
def define_problem(self):
y1,y2=var(["y1","y2"]) #bind the variables here
K=self.Kmatrix # shorthand
#Adding two single oscillators, but passing the coupling
eqs=SingleHarmonicOscillator("y1",K[0][0]*y1+K[0][1]*y2)
eqs+=SingleHarmonicOscillator("y2",K[1][0]*y1+K[1][1]*y2)
eqs+=InitialCondition(y1=1,y2=0) #Setting initial condition
eqs+=ODEFileOutput()
self.add_equations(eqs@"coupled_oscillator")
Method III
The final method is very similar to the second method, but the equations are not combined to a single equation on the domain "coupled_oscillator", but we use two different domains:
def define_problem(self):
#bind the variables. Both are called just "y", but on different domains
y1=var("y",domain="oscillator1") #note the important keyword argument domain!
y2=var("y",domain="oscillator2")
K=self.Kmatrix # shorthand
#Define the first harmonic oscillator
eqs1=SingleHarmonicOscillator("y",K[0][0]*y1+K[0][1]*y2)
eqs1+=InitialCondition(y=1)
eqs1+=ODEFileOutput()
#Define the second harmonic oscillator
eqs2=SingleHarmonicOscillator("y",K[1][0]*y1+K[1][1]*y2)
eqs2+=InitialCondition(y=0)
eqs2+=ODEFileOutput()
#Add both to different domains
self.add_equations(eqs1@"oscillator1"+eqs2@"oscillator2")
Note that the unknowns of both oscillators have the same name, namely "y". However, since we are creating two distinct equations on different domains ("oscillator1" and "oscillator2"), the same names do not interfere. However, we have to be careful, since the unknown obtained by var("y") will have different meanings in both domains, namely the unknown "y" defined in the particular domain. When coupling unknowns of different domains, it is therefore required to pass the keyword domain in the var() function to remove this ambiguity.
Also, the initial conditions are now separated, since initial conditions can only be set on the domain where the corresponding unknown is defined. Since there are also two instances of ODEFileOutput, each domain will write its own output file, containing just the single degree of freedom within this domain. The calculated solution is depicted in Fig. 3.3.
Fig. 3.3 Output for the user-defined harmonic oscillatior equation with damping and driving.