IndexError: index 1 is out of scope for axis 0 with size 1 / ForwardEuler

Question

I numerically solve for x (t) for a system of first order differential equations. System:

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

I applied the Forward Euler method to solve this problem as follows: Here is my code:

 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() 

I get the following Traceback:

 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 

I don’t understand what is probably wrong. I have already raised my eyes after resolving issues, but it does not help me. Can you find my mistake? I use the following code as a guide: 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 ) 

My goal is to build a schedule.

Answer 1

The problem, as Traceback says, comes from the line x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] ) . Let me replace it in my context:

• x is an array equal to [x0 * n], so its length is 1
• you repeat from 0 to n-2 (n doesn't matter here), and I am the index. In the beginning, everything is fine (there is no beginning here, apparently ... :(), but as soon as i + 1 >= len(x) <=> i >= 0 , the element x[i+1] does not exist. Here this the element doesn't exist since the beginning of the for loop.

To solve this problem, you should replace x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] ) with x.append(x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )) .

Answer 2

Problem with your line

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

Here you define x as a single element array of -200.0. You can do it:

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

or that:

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

Note: your import is rather confused. You import numpy modules three times into the header, and then import pylab (which already contains all numpy modules). If you want easy, with one single

 from pylab import * 

at the top you can use all the modules you need.