Bottom Page

• 0 Vote(s) - 0 Average
• 1
• 2
• 3
• 4
• 5
 Odeint to solve Mukhanov equation Messier087 Unladen Swallow Posts: 4 Threads: 2 Joined: Feb 2020 Reputation: 0 Likes received: 0 #1 Feb-12-2020, 02:12 PM 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, -(3-ep)*Y -((k**2)/H)*Y] ########################################## # 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  Posts: 673 Threads: 1 Joined: Mar 2018 Reputation: 90 Likes received: 103 #2 Feb-13-2020, 03:36 AM 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 Unladen Swallow Posts: 4 Threads: 2 Joined: Feb 2020 Reputation: 0 Likes received: 0 #3 Feb-16-2020, 11:02 PM 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  Posts: 673 Threads: 1 Joined: Mar 2018 Reputation: 90 Likes received: 103 #4 Feb-16-2020, 11:46 PM 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, -(3-ep_func(t))*Y -((k**2)/H_func(t))*Y]````H_func` and `ep_func` are arbitrary and can include some interpolations as well Messier087 Unladen Swallow Posts: 4 Threads: 2 Joined: Feb 2020 Reputation: 0 Likes received: 0 #5 Feb-17-2020, 05:05 PM Thank you so much! That was really helpful :D « Next Oldest | Next Newest »

Top Page

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

Forum Jump:

Users browsing this thread: 1 Guest(s)