如何在Python 3中使用龙格-库塔-4调用数值结果来集成ODE

Question

我正试图用数值方法解决(对于m_0)以下常微分方程:

dm0/dx=(((1-x)*(x*(2-x))**(1.5))/(k+x)**2)*(((x*(2-x))/3.0)*(dw/dx)**2 + ((8*(k+1))/(3*(k+x)))*w**2)

w和dw/dx的值已经用龙格-库塔四阶得到，k是一个固定的因子。我写了一段代码，从外部文件中调用w和dw/dx的值，然后把它们组织成一个数组，然后在函数中调用这个数组，然后运行集成。我的结果不像预期的那样，我不知道哪里出了问题。如果有人能帮我一把，我将不胜感激。谢谢你！

from math import sqrt
from numpy import array,zeros,loadtxt
from printSoln import *
from run_kut4 import *
m = 1.15                         # Just a constant.
k = 3.0*sqrt(1.0-(1.0/m))-1.0    # k in terms of m.
omegas = loadtxt("omega.txt",float)    # Import values of w
domegas = loadtxt("domega.txt",float)  # Import values of dw/dx
w = []    # Defines the array w to store the values w^2
s = 0.0
for s in omegas:
    w.append(s**2)       # Calculates the values w**2
omeg = array(w,float)    # Array to store the value of w**2
dw = []      # Defines the array dw to store the values dw**2
t = 0.0
for t in domegas:
    dw.append(t**2)    # Calculates the values for dw**2
domeg = array(dw,float)     # Array to store the values of dw**2
x = 1.0e-12              # Starting point of integration
xStop = (2.0 - k)/3.0    # Final point of integration 
def F(x,y):      # Define function to be integrated
    F = zeros(1)
    for i in domeg:  # Loop to call w^2, (dw/dx)^2       
        for j in omeg:
            F[0] = (((1.0-x)*(x*(2.0-x))**(1.5))/(k+x)**2)*((1.0/3.0)*x* (2.0-x)*domeg[i] + (8.0*(k+1.0)*omeg[j])/(3.0*(k+x))) 
            return F
y = array([((32.0*sqrt(2.0)*(k+1.0)*(x**2.5))/(15.0*(k**3)))])  # Initial condition for m_{0}
h = 1.0e-5          # Integration step
freq = 0            # Prints only initial and final values
X,Y = integrate(F,x,y,xStop,h)      # Calls Runge-Kutta 4
printSoln(X,Y,freq)                 # Prints solution

Answer 1

解释您的口头描述，有一个ODE用于omega, w'=F(x,w)，和一个耦合ODE用于m0, m'=G(x,m,w,w')。几乎总是解决这个问题的最佳方法是将其视为ODE系统，

def ODEfunc(x,y)
    w,m = y
    dw = F(x,w)
    dm = G(x,m,w,dw)
    return np.array([dw, dm])

，然后您可以将其插入到您选择的ODE求解器中，例如，虚构的

ODEintegrate(ODEfunc, xsamples, y0)