numpy 二维阵列中多个谐振子的编码四阶龙格库塔方法,卡在单个振荡器中[关闭]

wooyq4lh 于 11个月前发布在其他

关注(0)|答案(1)|浏览(85)

已关闭。此问题需要更多focused。目前不接受回答。
**要改进此问题吗？**更新此问题，使其仅针对editing this post的一个问题。

13天前关闭
Improve this question的
我之前在Python中用四阶龙格库塔方法数值求解了一般（非线性）pengulum。这里有一个与此相关的问题，Attempt to solve the nonlinear pendulum 2nd order differential equation using 4th order Runge-Kutta method, not getting expected result，它包含了代码和过程的细节。
我也试着用同样的方法解决一个由三个微分方程组成的“自旋”系统，也是成功的。
但是，这两个问题都涉及单个粒子。

**现在如果我有很多粒子。**所以我需要为每个粒子（没有相互作用）求解相同的微分方程。也许我可以尝试为所有粒子编写相同的步骤，但这很麻烦。而且，在这种情况下，我不能改变粒子数量，我必须重写整个过程。

我不是很熟悉 NumPy，我也迫不及待地想先详细了解它。
假设对于 * 单粒子/振荡器 *，我们有以下函数：

def f1(t,x,y): return y

def f2(t,x,y): return -k*sin(x)

字符串
然后是初始值：

k=1.0                      #parameter
t,x,y=0,8.0*pi/9.0,0       #initial values (t: second, x: radian, y: radian/second)
h=0.01                     #increment in t

型
RK 4循环：

T,X,Y=[t],[x],[y]          #lists to store data

# Loop:
for i in range(2000):
    a1=h*f1(t,x,y)                      
    b1=h*f2(t,x,y)                      
    a2=h*f1(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    b2=h*f2(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    a3=h*f1(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    b3=h*f2(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    a4=h*f1(t+h,x+a3,y+b3)              
    b4=h*f2(t+h,x+a3,y+b3)              
    x=x+(1/6)*(a1+2*a2+2*a3+a4)         # apprx. value of x1 for t+h
    y=y+(1/6)*(b1+2*b2+2*b3+b4)         # apprx. value of y1 for t+h
    t=t+h                               # current value of independent variable t
    T.append(t)
    X.append(x)
    Y.append(y)

型
现在，假设，* 我们在一个一维阵列中有两个粒子/振荡器，或者，假设，在一个二维阵列中有2x2=4个振荡器，该怎么办？*
我们如何才能通过使用基本的numpy技术来完成我们的任务？
我想到了如下的东西：

def f1(t,x[i][j],y[i][j]): return y[i][j]

def f2(t,x[i][j],y[i][j]): return -k*sin(x[i][j])

k=1.0                      #parameter
t,x[i][j],y[i][j]=0,8.0*pi/9.0,0       #initial values (t: second, x: radian, y: radian/second)
h=0.01                     #increment in t

型
但是有几个错误被显示出来，我知道我错过了很多东西。
这就是为什么我根本无法进行循环部分。
需要添加哪些内容？

numpy

来源：https://stackoverflow.com/questions/77710244/coding-4th-order-runge-kutta-method-for-multiple-harmonic-oscillators-in-a-2d-ar

1条答案

按热度按时间

vlurs2pr1#

与其处理多维数组，你可以使用一个大小为2N的一维numpy数组（其中N是振荡器的数量）。
在这里，您可以将位置存储在y[0], y[2], y[4], ...中，将速度存储在y[1], y[3], y[5], ...中
这样做的好处是可以保持完全矢量化的龙格-库塔例程。

import numpy as np
import matplotlib.pyplot as plt
import math

# 4th-order explicit Runge-Kutta
def rk4( x, y, dx, f ):
    dy1 = dx * f( x           , y             )
    dy2 = dx * f( x + 0.5 * dx, y + 0.5 * dy1 )
    dy3 = dx * f( x + 0.5 * dx, y + 0.5 * dy2 )
    dy4 = dx * f( x +       dx, y +       dy3 )
    return x + dx, y + ( dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4 ) / 6.0

# Equation parameters
N = 3                                            # number of systems
k = np.array ( [ 1.0, 4.0, 9.0 ] )               # stiffnesses
m = np.array ( [ 1.0, 1.0, 1.0 ] )               # masses
x0 = np.array( [ 0.5, 0.6, 0.7 ] )               # initial displacements
v0 = np.array( [ 0.0, 0.0, 0.0 ] )               # initial velocities

# Derivative function (outputs a numpy array)
def deriv( t, y ):
    f = np.zeros_like( y )
    f[0:2*N-1:2] = y[1:2*N: 2]                   # derivative of positions
    f[1:2*N  :2] = -k * y[0:2*N-1:2] / m         # "true" harmonic oscillator; use np.sin() otherwise
    return f

# Initialise
t = 0
y = np.array( [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ] )
y[0:2*N-1:2] = x0                                # set positions in elements 0, 2, 4, ...
y[1:2*N  :2] = v0                                # set velocities in elements 1, 3, 5, ...
nt = 400
dt = 0.01

# Initialise for plotting
tvals=[]
xvals=[[] for i in range(N)]                     # careful!!!

# Successive timesteps
tvals.append( t )
for i in range( N ): xvals[i].append( y[2*i] )
for i in range( nt + 1 ):
    t, y = rk4( t, y, dt, deriv )               # Runge-Kutta update
    tvals.append( t )
    for i in range( N ): xvals[i].append( y[2*i] )

# Plot
for i in range( N ): plt.plot( tvals, xvals[i] )
plt.show()

字符串
x1c 0d1x的数据

赞(0）回复(0）举报 11个月前

我来回答

numpy 二维阵列中多个谐振子的编码四阶龙格库塔方法,卡在单个振荡器中[关闭]

1条答案

相关问题

热门标签

最新问答