0

我正在尝试以数值方式求解（对于 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 的值已经使用 Runge-Kutta 4 阶在数值上找到了，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
4

1 回答 1

0

解释您的口头描述，对于 omega, 有一个 ODE，w'=F(x,w)对于 m0, 有一个耦合 ODE 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)
于 2015-09-08T14:13:07.800 回答