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 の値は、ルンゲクッタの 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,w'=F(x,w)の ODE と 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 に答える