python - Mukhaov 方程的解
问题描述
我目前正在尝试寻找幂律势的 Mukhanov 方程的解。我的代码目前看起来像这样,
#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
import math
plt.rcParams["font.size"] = 6
plt.autoscale(enable=True, axis='both', tight=None)
#Define Parameters.
n = float(input('Enter the value of n: '))
k_a = np.geomspace(start = 0.005, stop = 10, num = 1000)
t = np.linspace(start = 60.0, stop = 0.0, num = 1000) #t is the efold number.
k_0 = 0.005
k_max = 1e16
#########################################
#Slow-roll parameters (using phi^2 for now and t is the e-fold number)
def ep_func(t, n):
return (n**2)/(4*n*t + n**2)
def H_func(t, n):
return ((2*n*t + (n**2)/2)**(n/2))/3
def eta_func(t, n):
return 2*((2*n - n**2)/(2*n*t + ((n**2)/2)))
def a_func(t):
return np.exp(t)
##########################################
#n_s and r
def r_func(t, n):
return (16*(n**2))/(4*n*t + n**2)
def n_func(t, n):
return 1 -(6*(n**2)/(4*n*t + n**2)) + (2*(n*(n-1)))/(2*n*t + (n**2)/2)
##########################################
#Create the two simplified ODE's.
#This is in terms of the e-fold number not confromal time.
def sol(Y, t, k, ep, eta, H, a):
return np.array([Y[1], -(3 - ep_func(t, n) + eta_func(t, n))*Y[1] -((k**2)/((a_func(t)**2)*H_func(t, n)**2))*Y[0]])
##########################################
# Find the solution using Scipy odeint function while defining the I.C's.
Yr_0 = [1/(np.sqrt(2*k_max)), 0]
Yi_0 = [0, -np.sqrt(k_max/2)]
solutions = []
# Real part of the Mukhanov equation.
for k in k_a:
asolr = odeint(sol, Yr_0, t, args=(k, ep_func(t, n), eta_func(t, n), H_func(t, n), a_func(t)))
solutions.append(asolr)
# Imaginary Part of the Mukhanov equation.
for k in k_a:
asoli = odeint(sol, Yi_0, t, args=(k, ep_func(t, n), eta_func(t, n), H_func(t, n), a_func(t)))
solutions.append(asoli)
Yr = asolr[:,0]
Yi = asoli[:,0]
##########################################
#Plot the solution to the mukhanov equation.
def mukh(Yi):
return np.absolute(Yi)
m = mukh(Yi)
plt.plot(t, m, 'r-')
plt.yscale('log')
##########################################
#Define the the power spectrum.
#Mukhanov Power spectrum
def m_spec(k_a, Yr, Yi):
return 2.1e-9*((k_a**3)/0.005)*(np.absolute(Yr)**2 + np.absolute(Yi)**2)
#Slow-roll power spectrum
def s_spec(k_a, t, n):
return 2.1e-9*(k_a/0.005)**(-(2*n**2 + n)/(2*n*t + (n**2)/2))
pm = m_spec(k_a, Yr, Yi)
ps = s_spec(k_a, t, n)
#Plot the power spectrums and parameters in subplots
fig = plt.figure()
plt.subplot(2, 2, 1)
plt.plot(k_a, pm, 'r-')
plt.title('Power spectrum P(k) for n = ' +str(n))
plt.xlabel('wavenumber (k)')
plt.ylabel('P(K)')
plt.yscale('log')
plt.xscale('log')
plt.subplot(2, 2, 2)
plt.plot(k_a, ps, 'b-')
plt.title('Slow-roll Power spectrum P(k) for n = ' +str(n))
plt.xlabel('wavenumber (k)')
plt.ylabel('P(K)')
plt.yscale('log')
plt.xscale('log')
plt.subplot(2, 2, 3)
plt.plot(t, H_func(t, n), 'y-')
plt.title('Epsilon Vs e-fold')
plt.xlabel('e-fold')
plt.ylabel('$H$(N)')
#plt.yscale('log')
plt.subplot(2, 2, 4)
plt.plot(t, n_func(t, n), 'g-')
plt.title('$n_s$ Vs e-fold')
plt.xlabel('e-fold')
plt.ylabel('$n_s$')
#plt.yscale('log')
plt.tight_layout()
plt.show()
我已经仔细检查了微分方程及其初始条件是正确的。我相信问题在于解决方案的虚部。我目前得到一个看起来像这样的输出。
问题是左上角的图像。我希望一条直线具有大约 1 的近似恒定梯度,永远不会超过 10^-9,显然情况并非如此。
解决方案
推荐阅读
- javascript - 为什么fetchmock restore会给出无与伦比的get?
- elasticsearch - 多匹配查询分页
- javascript - 减少和求和对象条目数组
- javascript - 如何通过 html 或 javascript 在 chrome 中打开翻译对话框
- mysql - 我应该使用视图进行固定条件查询吗?
- spring-boot - 使用 gradle 的 OutOfMemoryError
- python - 在熊猫中连接两个数据帧时出现索引错误
- javascript - 调用数据库后如何通过href链接为每一行创建一个弹出窗口
- android-studio - Android Studio 始终使用 Proguard 和混淆调试构建 - 无法调试
- python - 替换 find_pos?