Решение вектора ODE в Python - PullRequest
Новая
       37

Решение вектора ODE в Python

0 голосов
/ 04 декабря 2018

Я хотел бы поблагодарить любого, кто ответит на этот вопрос заранее, так как это, вероятно, глупый вопрос и пустая трата времени.

Данный код принимает электромагнитные силы на частицу и предположительно строит ее траекторию,но я не могу понять, как использовать RK4 с векторами.Чувствую, что моя функция настроена неправильно. 

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv, det, eig
from numpy.random import randn

def rk4(f, x0, dt):
    tn = 0
    xn = x0

    while True:        
        yield tn,xn

        k1 = dt*f(tn,xn)
        k2 = dt*f(tn+dt/2,xn+k1/2)
        k3 = dt*f(tn+dt/2,xn+k2/2)
        k4 = dt*f(tn+dt,xn+k3)

        xn = xn + (k1+2*k2+2*k3+k4)/6
        tn = tn + dt
#--------------------------------------------------------

def f(t,X):
    x,y,z,xv,yv,zv = X
    v = [xv,yv,zv]
    E = [k*x , k*y , -2*z]
    a = (q/m)*(E+np.cross(v,B))
    Xdot = np.array([v,a])
    return Xdot

q , m , B , k = 1.6e-19, 40, [0,0,1], 10e+4


X0 = np.array([0.001 , 0 , 0.001 , 100 , 0 , 0]) 
solver = rk4(f,X0,10e-7) 

ts = []
Xs = []
for t,X in solver:
    ts.append(t)
    Xs.append(X[0])
    if t > 0.001:
        break


#Xs = np.array(Xs) 
#plt.figure()
#plt.plot(ts,Xs)

Буду очень признателен за любые советы или подсказки. Я подозреваю, что проблема связана с тем, как я настроил функцию, так как код вылетает, когда он пытается его выполнить.RK4.

1 Ответ

0 голосов
/ 04 декабря 2018

Вы не можете применять векторную арифметику к простым спискам Python, поэтому сначала преобразуйте списки в массивы.Вам нужно вернуть плоский вектор, а не матрицу, поэтому используйте объединение массивов для объединения двух частей. 

def f(t,X):
    x,y,z,xv,yv,zv = X
    v = np.array([xv,yv,zv])  #  or v = X[3:]
    E = np.array([k*x , k*y , -2*z]) # or E=k*X[:3]; E[2]=-2*X[2]
    a = (q/m)*(E+np.cross(v,B))
    Xdot = np.concatenate([v,a])
    return Xdot

Последнее изменение см. Объединение двух одномерных массивов NumPy

...