IndexError: 索引1超出了大小为1的轴0的范围/向前欧拉法

Question

我正在求解一组一阶微分方程的x(t)数值解，该系统为：

dy/dt=(C)\*[(-K\*x)+M*A]

我已经按照以下方式实现了向前欧拉法来解决这个问题。这是我的代码：

import matplotlib
import numpy as np
from numpy import *
from numpy import linspace
from matplotlib import pyplot as plt


C=3
K=5
M=2
A=5
#------------------------------------------------------------------------------
def euler (f,x0,t):
    n=len (t)
    x=np.array ([x0*n])
    for i in xrange (n-1):
        x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )
    return x



#---------------------------------------------------------------------------------          
if __name__=="__main__":
    from pylab import *

    def f(x,t): 
        return (C)*[(-K*x)+M*A]

    a,b=(0.0,10.0)
    n=200
    x0=-1.0
    t=linspace (a,b,n)

    #numerical solutions
    x_euler=euler(f,x0,t)

    #compute true solution values in equal spaced and unequally spaced cases
    x=-C*K
    #figure
    plt.plot (t,x_euler, "b")
    xlabel ()
    ylabel ()
    legend ("Euler")

    show()
`
M=2
A=5
#----------------------------------------------------------------------------
def euler (f,x0,t):
    n=len (t)
    x=np.array ([x0*n])
    for i in xrange (n-1):
        x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )
    return x



#---------------------------------------------------------------------------          
if __name__=="__main__":
    from pylab import *

    def f(x,t): 
        return (C)*[(-K*x)+M*A]

    a,b=(0.0,10.0)
    n=200
    x0=-1.0
    t=linspace (a,b,n)

    #numerical solutions
    x_euler=euler(f,x0,t)

    #compute true solution values in equal spaced and unequally spaced cases
    x=-C*K
    #figure
    plt.plot (t,x_euler, "b")
    xlabel ()
    ylabel ()
    legend ("Euler")

    show()

我遇到了以下错误信息：

Traceback (most recent call last):
  File "C:/Python27/testeuler.py", line 50, in <module>
    x_euler=euler(f,x0,t)
  File "C:/Python27/testeuler.py", line 28, in euler
    x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )
IndexError: index 1 is out of bounds for axis 0 with size 1

我不明白可能出了什么问题。我已经查找了解决方案，但这并没有帮助我。你能找到我的错误吗？ 我使用以下代码作为参考： def euler( f, x0, t ):

    n = len( t )
    x = numpy.array( [x0] * n )
    for i in xrange( n - 1 ):
        x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )

    return x
if __name__ == "__main__":
    from pylab import *

    def f( x, t ):
        return x * numpy.sin( t )

    a, b = ( 0.0, 10.0 )
    x0 = -1.0

    n = 51
    t = numpy.linspace( a, b, n )

    x_euler = euler( f, x0, t )

我的目标是绘制这个函数。

- Mlle Blanche

---

x=np.array([x0*n])产生一个只有一个元素的数组。你当时想要做什么？ - mdurant

Answer 1

问题就像Traceback所说的那样，来自于这一行代码：x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )。让我们将它放回到上下文中：

• x是一个等于[x0 * n]的数组，因此它的长度为1
• 你正在从0到n-2进行迭代（这里n并不重要），i是索引。在开始时，一切都很好（这里似乎没有开始... :( ），但是一旦i + 1 >= len(x) < => i >= 0，元素x[i+1]不存在。在这里，自for循环开始以来，这个元素就不存在。

要解决这个问题，您必须将x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )替换为x.append(x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] ))。

- Spirine

Answer 2

问题出在你的线路上

x=np.array ([x0*n])

在这里，你将x定义为一个-200.0的单元素数组。你可以这样做：

x=np.array ([x0,]*n)

或者这个：

x=np.zeros((n,)) + x0

注意：您的导入相当混乱。您在标题中导入了numpy模块三次，然后稍后导入pylab（它已经包含所有numpy模块）。如果您想简单一些，只需使用一个单一的

from pylab import *

在顶部的行中，您可以使用所有需要的模块。

- leeladam