Odeint to solve Mukhanov equation

[Messier087]

Hi,

I am currently trying to solve the Mukanov equation which contains variables that change with time. This requires me to loop the odeint function but I'm getting errors. The code is below.

#Code to compute the power spectrum from the solution to the Mukhanov equation.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

plt.rcParams["font.size"] = 8
plt.autoscale(enable=True, axis='both', tight=None)

#Define Parameters. 
k = 10
t = np.linspace(0.0, 60.0, 10000) #t is the efold number.

##########################################

#Slow-roll parameters (phi^2 for now)

def epsilon(t):
    return 1/(2*t + 1)

def H(t):
    return ((4*t + 2)**2)/3

ep = epsilon(t),

H = H(t)

##########################################
#Create the two simplified ODE's.
#This is in terms of the e-fold number not confromal time. 
def sol(Y, t, k, ep, H):
    return [Y[1], -(3-ep)*Y[1] -((k**2)/H)*Y[0]]
##########################################

# Find the solution using Scipy odeint function while defining the I.C's. 
Y0 = [1e-3, 0]

solutions = []
for ep, H in range(len(t)):    
    asol = odeint(sol, Y0, t, args=(k, ep, H), h0 = 0.005)  
    solutions.append(asol)

y = asol[:,0]
##########################################

#Plot the solution to the mukhanov equation. 
#plt.plot(t, y, 'r-')
#plt.yscale('log')

##########################################

#Define the the power spectrum.
def m_spec(k, y):
    return ((k**3)/2*np.pi**2)*np.absolute(y)**2

pm = m_spec(k, y)

#Plot the power spectrum.
for asol in solutions:
    plt.plot(t, pm, 'r-')

plt.yscale('log')   
plt.show()

The problem is in the H term. If I explicitly but H into the solution function and just loop over ep it works, but I cannot loop over both and I'm not sure why. I'm really new to python so I could be doing something stupid!

[scidam]

As far as I understood, the problem is in lines 38-41. You can resolve it, e.g. assuming that ep = H,

solutions = []
for ep, H in zip(range(len(t)), range(len(t))):    
    asol = odeint(sol, Y0, t, args=(k, ep, H), h0 = 0.005)  
    solutions.append(asol)

Another approach would be to use dictionary to store results:

solutions = dict()
ep_values = [1,2,3]
H_values =[4,3,1]
for ep in ep_values:
    for H in H_values:
        asol = odeint(sol, Y0, t, args=(k, ep, H), h0 = 0.005)  
        solutions.update({(ep, H): asol})

But this will require to restructure the rest part of your code.

[Messier087]

I've been doing some reading, is there any way of using interpolation functions to feed H and eq into odeint? They change over time and odeint doesn't appear to like that.

[scidam]

If H and ep depend on t you always can do the following:

def H_func(t):
    return ... 

def ep_func(t):
    return ... 

def sol(Y, t, k):
    return [Y[1], -(3-ep_func(t))*Y[1] -((k**2)/H_func(t))*Y[0]]

H_func and ep_func are arbitrary and can include some interpolations as well

[Messier087]

Thank you so much! That was really helpful :D