Решения Lane-Emden с использованием RK4 (жестко) - PullRequest
Новая
       18

Решения Lane-Emden с использованием RK4 (жестко)

0 голосов
/ 14 сентября 2018

[LE Solutions ] [1]

Я работаю над поиском решения для построения уравнения Лейна-Эмдена для значений n = [0,6] с интервалами1/2.Я новичок в Python, и не могу понять, как использовать RK4 для этой работы.Пожалуйста помоги!

Текущий прогресс.«TypeError: неподдерживаемые типы операндов для Pow: 'int' и 'list" в строке 37 в main.py "Ошибка появилась только после того, как я добавил в уравнения, определенные как r2, r3, r4 и k2, k3, k4. 

import numpy as np
import matplotlib.pyplot as plt

n = [0,1,2,3,4,5,6,7,8,9,10,11,12,13]
theta0 = 1
phi0 = 0
step = 0.01
xi0 = 0
xi_max = 100

theta = theta0
phi = phi0
xi = xi0 + step

Theta = [[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
Phi = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
Xi = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]

for i in n:
    Theta[i].append(theta)
    Phi[i].append(phi)
    Xi[i].append(xi)
 
def dTheta_dXi(phi,xi): #r1
    return -phi/xi**2
    
def r2(phi,xi):
    return dTheta_dXi(phi+step,xi+step*dTheta_dXi(phi,xi))

def r3(phi,xi):
    return dTheta_dXi(phi+step,xi+step*r2(phi,xi))
    
def r4(phi,xi):
    return dTheta_dXi(phi+step,xi+step*r3(phi,xi))

def dPhi_dXi(theta,xi,n): #k1
    return theta**(n)*xi**2
    
def k2(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*dPhi_dXi(theta,xi,n),n)

def k3(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*k2(theta,xi,n),n)
    
def k4(theta,xi,n):
    return dPhi_dXi(theta+step,xi+step*k3(theta,xi,n),n)
    
for i in n:	
    while xi < xi_max:
        if theta < 0:
            break
        dTheta = (step/6)*(dTheta_dXi(phi,xi)+2*r2(phi,xi)+2*r3(phi,xi)+r4(phi,xi))
        dPhi = (step/6)*(dPhi_dXi(theta,xi,i/2.)+2*k2(theta,xi,n)+2*k3(theta,xi,n)+k4(theta,xi,n))
        theta += dTheta
        phi += dPhi
        xi += step
        Theta[i].append(theta)
        Phi[i].append(phi)
        Xi[i].append(xi)
    print i/2., round(xi,4), round(dTheta_dXi(phi,xi),4), round(xi/3./dTheta_dXi(phi,xi),4), round(1./(4*np.pi*(i/2.+1))/dTheta_dXi(phi,xi)**2,4)
    theta = theta0
    phi = phi0
    xi = xi0 + step

1 Ответ

0 голосов
/ 15 сентября 2018

Вы должны еще раз изучить, какое обновление идет куда. xi является независимой переменной и, следовательно, получает только обновления 0.5*step и step, обновления theta используют производные dTheta_dXi и аналогично phi обновляется с использованием наклонов dPhi_dXi

def r2(phi,xi):
    return dTheta_dXi(phi+0.5*step*dPhi_dXi(theta,xi,n),xi+0.5*step)

def k2(theta,xi,n):
    return dPhi_dXi(theta+0.5*step*dTheta_dXi(phi,xi),xi+0.5*step,n)

def r3(phi,xi):
    return dTheta_dXi(phi+0.5*step*k2(theta,xi,n),xi+0.5*step)

и т.д..

Теперь можно видеть, что из-за связанной природы уравнения вам нужны и theta, и phi в качестве аргументов повсюду. Кроме того, даже если это сработает в конце, вы в конечном итоге вычисляете много значений несколько раз, где для сборки всего в одном цикле требуется только одно вычисление.

...