2016-10-10 266 views
0

我想解决这个方程:解决使用odeint二阶ODE在Python

Y '' + A Y” - B y = 0的

其中y,A和B是的功能同样的变量“a”

我尝试下面的代码:

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


x0 = 0.,0.1 # initial conditions 
oe = 0.0 # cosmological constant density parameter 
om = 1. # matter density parameter 
h = 2.3E-18 # Hubble constant 
w = -1. 

a = np.linspace(0.1,1,10) # scale factor 


def H(a): # Hubble rate equation 
    return h*np.sqrt((om/(a**3.))+(oe/(a**(3.*(1.+w))))) 

def A(a): # differential equation term 
    return (-3./(2.*a)) 

def B(a): # differential equation term 
    return (3.*om*(h**2.))/((2.*(a**5.))*(H(a)**2.)) 



def system(X,a): # differential equation system 
    X0 = X[0] 
    X1 = X[1] 
    X2 = -A(a)*X1 + B(a)*X0 

    return X1,X2 



x = odeint(system,x0,a) 



pyplot.semilogx(a,x, linestyle='-', c="k", linewidth="2") 

它返回一个情节是没有意义的。我应该只得到一个情节,最大值“x”在a = 1时为1。但我得到了以下情节:

情节,我得到:

Plot that I get

和预期的结果是一样的连续线下图中:

enter image description here

任何建议?

+0

那你得到了什么?为什么不添加一个数字?请阅读[问]。另外,H^2对数值精度不会有好处。如果可能,我建议重写您的方程,使得数量在1的数量级上, –

+0

只需按照示例一步一步http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate。 odeint.html – outoftime

+0

该代码的工作原理应该如此,您应该将问题转移到物理论坛,讨论要更改的模型以使模型符合您的期望。 – LutzL

回答

0

请尝试在功能系统(X,a)中的B(a)* X0之前更改符号。根据你的初始方程,它应该是负值。

+0

我的不好。原方程是y“+ A * y' - B * y –