在 Python 中求解向量 ODE

Question

我要感谢提前回答这个问题的人，因为这可能是一个愚蠢的问题，回答是浪费时间。

给定的代码采用粒子上的电磁力并据称绘制了它的轨迹，但我不知道如何将 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)

我非常感谢任何提示或提示我怀疑问题源于 iv 设置函数的方式，因为一旦代码尝试通过 rk4 运行它就会崩溃。

Answer 1

您不能将向量运算应用于简单的 python 列表，因此请先将列表转换为 numpy 数组。您需要 return 一个平面向量，而不是矩阵，因此使用数组连接来连接两个部分。

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

最后的变化见Concatenating two one-dimensional NumPy arrays

在 Python 中求解向量 ODE

Solving Vector ODE in Python

python

physics

vector

ode

python-3.x