Bottom Page

Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
 Odeint to solve Mukhanov equation
#1
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! Shifty
Quote
#2
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.
Quote
#3
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.
Quote
#4
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
Quote
#5
Thank you so much! That was really helpful :D
Quote

Top Page

Possibly Related Threads...
Thread Author Replies Views Last Post
  Mukhanov equation + odeint Messier087 0 355 Mar-28-2020, 04:03 PM
Last Post: Messier087
  Can't find a way to solve this nonlinear equation Alex009988 2 987 Aug-16-2019, 01:50 AM
Last Post: scidam
  odeint to solve Schrodinger equation kiyoshi7 14 5,072 Nov-23-2018, 11:49 AM
Last Post: kiyoshi7
  python odeint keeps giving me size of array error kiyoshi7 1 1,586 Nov-01-2018, 02:03 AM
Last Post: j.crater

Forum Jump:


Users browsing this thread: 1 Guest(s)