Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #simple exponential potential
- # u = K*phi'/H0; v = omega_matter**(1/3)*(1+z); w = l*K*phi' - ln((K**2)*V0/H0**2)
- # f,g,h are functions of derivation of u,v,w respectively derieved w.r.t t*H0 = T
- import matplotlib.pyplot as plt
- import numpy as np
- import math
- def f(T,u,v,w):
- return -3*u*((v**3 + (u**2)/6 + np.exp(-w)/3)**0.5) + l*np.exp(-w)
- def g(T,u,v,w):
- return -v*(v**3 + (u**2)/6 + np.exp(-w)/3)**0.5
- def h(T,u):
- return l*u
- p = 0.1
- q = 1.0
- dh = 0.1
- n = (q-p)/dh
- u = [0.0]
- v = [1100]
- T = [0.00001]
- w = [-44.927]
- l = 1.3
- for i in range(0,int(n)):
- k1 = f(T[i],u[i],v[i],w[i])
- r1 = g(T[i],u[i],v[i],w[i])
- print k1, r1
- s1 = h(T[i],u[i])
- print s1
- k2 = f(T[i] + 0.5*dh,u[i] + k1*dh,v[i] + k1*dh,w[i] + k1*dh)
- r2 = g(T[i] + 0.5*dh,u[i] + r1*dh,v[i] + r1*dh,w[i] + r1*dh)
- s2 = h(T[i] + 0.5*dh,u[i] + s1*dh)
- print k2,r2,s2
- k3 = f(T[i] + 0.5*dh,u[i] + k2*dh,v[i] + k2*dh,w[i] + k2*dh)
- r3 = g(T[i] + 0.5*dh,u[i] + r2*dh,v[i] + r2*dh,w[i] + r2*dh)
- s3 = h(T[i] + 0.5*dh,u[i] + s2*dh)
- k4 = f(T[i] + dh,u[i] + dh*k3,v[i] + dh*k3,w[i] + k3*dh)
- r4 = g(T[i] + dh,u[i] + r3*dh,v[i] + dh*r3,w[i] + r3*dh)
- s4 = h(T[i] + dh,u[i] + dh*s3)
- T == T.append(T[i] + dh)
- u == u.append(u[i] + (dh/6)*(k1 + 2.0*k2 + 2.0*k3 + k4))
- v == v.append(v[i] + (dh/6)*(r1 + 2.0*r2 + 2.0*r3 + r4))
- w == w.append(w[i] + (dh/6)*(s1 + 2.0*s2 + 2.0*s3 + s4))
- plt.plot(T,u, '-b')
- plt.plot(T,v, '-r')
- plt.plot(T,w, '-g')
- plt.title('quintessence cosmological model')
- plt.show()
Add Comment
Please, Sign In to add comment