2016-05-24 22 views
3

Ich versuche, ein Feder-Masse-Problem sowohl mit der Euler- als auch der Range-Kutta-Methode zu lösen und die Kurven zu vergleichen. Ich habe Funktionen für Euler und Runge-Kutta geschrieben, aber nachdem ich die Funktion für mein Problem aufgerufen habe, scheint es, dass meine Handlung keine Daten zeigt. Bitte helfen Sie mir die Handlung zu beheben und prüfen, ob es ein Fehler in meinem Code ist dankVersuch, DE zweiter Ordnung durch Euler und Runge_Kutta-Methode zu lösen

#function Euler 
def euler (y, t, dt, derivative): 
    y_next = y + derivative(y, t) * dt 
    return y_next 

# function Runge-Kutta 
# 2nd order Runge-Kutta method routine 

def Runge_Kutta (y, time, dt, derivative): 
    k0 = dt * derivative (y, time) 
    k1 = dt * derivative (y + k0, time + dt) 
    y_next = y + 0.5 * (k0 + k1) 
    return y_next 

und hier ist das Problem, das ich

[![""" A spring and mass system. the coefficient of friction \mu is not negligible.generate a position vs. time plot for the motion of the mass, given an initial displacement x = 0.2m , spring constant k = 42 N/m , mass m =0.25 Kg, coefficient of friction \mu = 0.15 and initial velocity v = 0 

F = -kx +/-mu mg """ 


from pylab import * 
from Runge_Kutta_routine import Runge_Kutta 
from eulerODE import euler 

N = 500  #input ("How many number of steps to take?") 
x0 = 0.2 
v0 = 0.0 
tau = 3.0  #input ("What is the total time of the simulation in seconds?") 
dt = tau /float (N-1) 

k = 41.0  #input (" what is the spring constant?") 
m = 0.25  #input ("what is the mass of the bob?") 
gravity = 9.8 
mu = 0.15  #input ("what is the coefficient of friction?") 

""" we create a Nx2 array for storing the results of our calculations. Each 2- element row will be used for the state of the system at one instant, and each instant is separated by time dt. the first element in each row will denote position, the second would be velocity""" 

y = zeros (\[N,2\]) 

y \[0,0\] = x0 
y \[0,1\] = v0 

def SpringMass (state, time): 
    """ we break this second order DE into two first order DE introducing dx/ dt = v & dv/dt = kx/ m +/- mu g.... 

Note that the direction of the frictional force changes depending on the sign of the velocity, we handle this with an if statement.""" 


    g0 = state\[1\] 
    if g0 > 0: 
     g1 = -k/m * state \[0\] - gravity * mu 
    else: 
     g1 = -k/m * state \[0\] + gravity * mu 

    return array (\[g0, g1\]) 

# Now we do the calculations 
# loop only N-1 so that we don;t run into a problem addresssing y\[N+1\] on the last point 

    for j in range (N-1): 
     #y \[j+1\] = euler (y\[j\] , 0, dt, SpringMass) 
     y \[j+1\] = Runge_Kutta (y\[j\], 0 , dt, SpringMass) 

# Now we plot the result 

time = linspace (0 , tau, N) 
plot (time, y\[:,0\], 'b-', label ='position') 

xlabel('time') 
ylabel('position') 


show()][1]][1] 
+0

Haben Sie versucht, die Werte zu drucken, z. um zu sehen, dass sie nicht NaN oder Unendlichkeit sind? – user20160

Antwort

1

Es sieht aus wie Ihre Schleife mit for j in range (N-1): beginnen zu lösen versuchen das berechnet Array y ist eingerückt, so betrachtet Python diese Zeilen als Teil der Funktion SpringMass. Da diese Zeilen nach der return-Anweisung stehen, werden sie nie ausgeführt.

Um dies zu korrigieren, verschieben Sie diese Zeilen so, dass die Zeile for keine Einrückung hat und die anderen Zeilen nur vier Leerzeichen Einrückung haben. Sehen Sie, ob das Ihr Problem löst.

Beachten Sie, dass Ihr hier beschriebener Code immer noch nicht funktioniert. Sie haben externe Backslashes vor eckigen Klammern, schreiben Sie Ihre Euler und Runge_Kutta Funktionen in einem einzelnen unbenannten Modul, aber der Hauptcode erwartet, dass sie in zwei verschiedenen Modulen usw. sind. Sie haben auch viele Beispiele für schlechten Stil. Diese Probleme sind wahrscheinlich, warum Sie noch keine (andere) Antwort bekommen haben. Tun Sie sich selbst einen Gefallen und reinigen Sie Ihren Code vor der Veröffentlichung hier und clean up your style.