1

我是 Python 的新手，我对编程语言的了解还处于初级阶段，所以我复制了此处显示的 Runge-Kutta Python 脚本并根据我的目的对其进行了修改。这是我当前的脚本：

import numpy as np
from matplotlib import pyplot as plt
a=0
b=np.pi
g=9.8
l=1
N=1000

def RK4(f):
    return lambda t, y, dt: (
            lambda dy1: (
            lambda dy2: (
            lambda dy3: (
            lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
            )( dt * f( t + dt  , y + dy3   ) )
        )( dt * f( t + dt/2, y + dy2/2 ) )
        )( dt * f( t + dt/2, y + dy1/2 ) )
        )( dt * f( t       , y         ) )

from math import sqrt
dy = RK4(lambda t, y: -y)

t, y, dt = 0., 1., np.divide((b-a),float(N))
i=0
T=np.zeros((N+1,1))
DY=T
Y=T
while t < (b-dt):
    T[i]=t
    DY[i]=dy(t,y,dt)
    t, y = t + dt, y + dy( t, y, dt )
    Y[i]=y
    i=i+1

plt.figure(1)
plt.plot(T,Y)
plt.show()

你可以忽略前几行中的 g 和 l 变量，我本来打算求解单摆的问题，但后来我记得这是一个一阶 ODE 求解器，所以现在我的 ODE 是dy/dx=-y. 我一直在 IPython 中运行它。我期待T成为一个数组，基本上相当于linspace(0,pi,N+1)MATLAB 中的数组。所以 T 是一组 N+1 均匀间隔的值，介于（包括）0 和 pi 之间，但这是其内容的样本（作为 running 的输出T）：

In [101]: T
Out[101]:
array([[ 0.99686334],
       [ 0.99373651],
       [ 0.9906195 ],
       ...,
       [ 0.04334989],
       [ 0.04321392],
       [ 0.        ]])

（包括输入和输出行，以提供一些关于我在这里指的是什么的上下文，以防不清楚）。哦，顺便说一句，如果你想知道为什么我没有使用T=np.linspace(a,b,num=N+1)它而不是在这个循环中定义它，那是因为这给出了类似不寻常的 T 数组。

4

1 回答 1

3

和

DY=T
Y=T

您只是在复制引用，它们都指向同一个数组对象。本质上，所有内容都将Y是最后分配的内容。

用于T.copy()获取单独的数组对象。

于 2016-01-01T00:44:11.047 回答