2015-05-05 13 views
6

Sto risolvendo numericamente per x (t) per un sistema di equazioni differenziali del primo ordine. Il sistema è:IndexError: index 1 è fuori limite per l'asse 0 con dimensione 1/ForwardEuler

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

Ho implementato il metodo Forward Eulero per risolvere questo problema nel modo seguente: Qui è il mio codice:

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

ottengo seguente 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 

Non capisco cosa sia probabilmente sbagliato. Ho già alzato gli occhi dopo aver risolto le domande, ma non mi aiuta. Scoprirai errore mio? Sto usando seguente codice come un orientamento: 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) 

Il mio obiettivo è quello di tracciare la funzione.

+1

'x = np.array ([x0 * n])' produce una matrice con un elemento. Cosa stavi cercando di fare? – mdurant

risposta

5

Il problema, come dice Traceback, deriva dalla linea x[i+1] = x[i] + (t[i+1] - t[i]) * f(x[i], t[i]). Diamo sostituirlo nel suo contesto:

  • x è un array uguale a [x0 * n], quindi la sua lunghezza è di 1
  • si sta scorrendo da 0 a n-2 (n non importa qui), e io sono l'indice. All'inizio, tutto è ok (qui non c'è inizio apparentemente ... :(), ma non appena i + 1 >= len(x) < =>i >= 0, l'elemento x[i+1] non esiste. Qui, questo elemento non esiste dall'inizio di . il ciclo for

per risolvere questo, è necessario sostituire x[i+1] = x[i] + (t[i+1] - t[i]) * f(x[i], t[i]) da x.append(x[i] + (t[i+1] - t[i]) * f(x[i], t[i]))

9

il problema è con la vostra linea

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

Qui si definisce x come una matrice singola voce di. - 200.0 potrebbe fare questo:

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

o questo:

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

Nota: le vostre importazioni sono abbastanza confusi. Importa moduli numpy per tre volte nell'intestazione e successivamente importa pylab (che contiene già tutti i moduli numpy). Se vuoi andare facile, con una singola

from pylab import * 

linea nella parte superiore è possibile utilizzare tutti i moduli necessari.

Problemi correlati