Simulation von Bewegungen

Häufig steht man in der Physik vor dem Problem, dass man zwar mithilfe des dritten Newtonschen Gesetzes eine Gleichung formulieren kann, die die Bewegung einer Masse in einem Kraftfeld genau beschreibt, doch ist diese Gleichung mit mathematischen Mitteln nicht lösbar.In diesem Fall greift man zu einer numerischen Methode, die zwar nur eine Approximation der wahren Lösung ergibt, doch bei genügendem Rechenaufwand genau genug für reale Fragestellungen ist.

Wir beginnen mit der Berechnung der Bahn, die ein abgeworfener Ball (eine Kugel) beschreibt. Zuerst der freie Fall im Vakuum, dann mit Luftreibung, schließlich den schiefen Wurf.

Als Lösungsmethode wählen wir ein Runge-Kutta-Verfahren mit einer festen Schrittweite von 1/100 Sekunde, was für unsere Ansprüche sicher genügt.

1. Freier Fall

## freier Fall ohne Reibung

from math import *


###########################################################
# Solution of y''(t) = f(t,y,y'), y(t0)=y0, y'(t0)=y'0
# RKN4 y''=f(t,y,y')

def rkn4_1(t,y,yd,f,h):
    k1 = f(t,y,yd)
    k2 = f(t+h/2.0,y+h/2.0*yd+h*h/8.0*k1,yd+h/2.0*k1)
    k3 = f(t+h/2.0,y+h/2.0*yd+h*h/8.0*k1,yd+h/2.0*k2)
    k4 = f(t+h,y+h*yd+h*h/2.0*k3,yd+h*k3)
    ynew = y+h*yd+h*h/6.0*(k1+k2+k3)
    ydnew = yd+h/6.0*(k1+2.0*k2+2.0*k3+k4)
    return (t+h,ynew,ydnew)


###########################################################
# freier Fall im Vakuum
# F = ma = -mg bzw. xd' = -g = f(t,x,xd) 

def f(t,x,v):
    ableitung = -g 
    return ableitung


        
def fall(hoehe=1000,geschwindigkeit=0):

    global g                    # für Funktion f verfügbar machen
    g = 9.81                    # Fallbeschleunigung

    t = 0.0                     # laufende Zeit, Startwert
    x = float(hoehe)            # Startwert x
    xd = float(geschwindigkeit) # Startwert xd
    h = 0.01                    # Integrations-Zeitschrittweite

    ## für hübsch formatierte Ausgabe
    printformat = "t = %4.2f, x = %8.4f, v = %8.4f"

    print printformat%(t,x,xd)  # Startsituation
    
    while x>=0.0:                   # solange über Erdoberfläche
        for i in range(100):        # alle 100 Schritte ausgeben
            t,x,xd = rkn4_1(t,x,xd, f,h)
            if x<=0.0: break 
        print printformat%(t,x,xd)

ein Testlauf für def Fall aus einer Höhe von 1000 Meter mit Anfangsgeschwindigkeit 0 liefert:

>>> fall(1000,0)

t = 0.00,  x = 1000.0000, v = 0.0000
t = 1.00,  x = 995.0950, v = -9.8100
t = 2.00,  x = 980.3800, v = -19.6200
t = 3.00,  x = 955.8550, v = -29.4300
t = 4.00,  x = 921.5200, v = -39.2400
t = 5.00,  x = 877.3750, v = -49.0500
t = 6.00,  x = 823.4200, v = -58.8600
t = 7.00,  x = 759.6550, v = -68.6700
t = 8.00,  x = 686.0800, v = -78.4800
t = 9.00,  x = 602.6950, v = -88.2900
t = 10.00, x = 509.5000, v = -98.1000
t = 11.00, x = 406.4950, v = -107.9100
t = 12.00, x = 293.6800, v = -117.7200
t = 13.00, x = 171.0550, v = -127.5300
t = 14.00, x =  38.6200, v = -137.3400
t = 14.28, x =  -0.2198, v = -140.0868
                            

Wir sehen das lineare Wachsen der Geschwindigkeit und die quadratisch wachsende Flugstrecke. Da wir immer 100 Rechenschritte mit 1/100 Sekunde zwischen den Ausgabezeilen durchführen, erhalten wir sekundenweise Zwischenergebnisse. Diese Werte könnte man allerdings, wie Du aus dem Unterricht weißt, auch leicht mathematisch mit Zettel, Bleistift und Integralrechnung erhalten. Aber jetzt kommt etwas, was analytisch nicht mehr lösbar ist:

2. Fall mit Luftreibung

Wir verwenden die Newton'sche Näherung einer Luftreibung, bei der die bremsende Kraft dem Quadrat der Geschwindigkeit des Körpers proportional ist. Nicht berücksichtigt wird dabei die Änderung der Reynoldszahl bei Geschwindigkeitsänderungen und die Veränderung der Situation bei Überschallgeschwindigkeit. Das würde uns zu weit führen. Außerdem nehmen wir keine Rücksicht auf die Oberflächenbeschaffenheit des Körpers (Golfball, Fußball,...) und eine eventuelle Drehbewegung (Spin) des Körpers.

Das Programm sieht so aus:

# Fall mit Luftreibung

from math import *


###########################################################
# Solution of y''(t) = f(t,y,y'), y(t0)=y0, y'(t0)=y'0
# RKN4 y''=f(t,y,y')

def rkn4_1(t,y,yd,f,h):
    k1 = f(t,y,yd)
    k2 = f(t+h/2.0,y+h/2.0*yd+h*h/8.0*k1,yd+h/2.0*k1)
    k3 = f(t+h/2.0,y+h/2.0*yd+h*h/8.0*k1,yd+h/2.0*k2)
    k4 = f(t+h,y+h*yd+h*h/2.0*k3,yd+h*k3)
    ynew = y+h*yd+h*h/6.0*(k1+k2+k3)
    ydnew = yd+h/6.0*(k1+2.0*k2+2.0*k3+k4)
    return (t+h,ynew,ydnew)


###########################################################
# Fall einer Kugel mit Newton'scher Reibungskraft
# F = ma bzw.  xd' = -g + cw.rho*A*v**2/m


def f(t,x,v):
    # die Kraft wirkt immer bremsend
    if v>=0.0: abl = -g -factor*v**2  
    else:      abl = -g +factor*v**2 
    return abl
        

def fall(hoehe=1000,geschwindigkeit=0,cw=0.43):

    cw = cw         # Widerstandsbeiwert
    m = 0.8         # Masse des Körpers
    r = 0.03        # ungefährer Radius des Körpers
    A = 4.0*pi*r**2 # Querschnittsfläche des Körpers
    rho = 1.213     # Luftdichte
    global g        # für f verfügbar machen
    g = 9.81        # Fallbeschleunigung
    global factor   
    factor = cw*rho*A/m

    t = 0.0                     # laufende Zeit
    x = float(hoehe)            # Startwert x
    xd = float(geschwindigkeit) # Startwert xd
    h = 0.01                    # Integration-Zeitschrittweite

    printformat = "t = %4.2f, x = %8.4f, v = %8.4f"
    print printformat%(t,x,xd)
    
    while x>=0.0:
        for i in range(100):
            t,x,xd = rkn4_1(t,x,xd, f,h)
            if x<=0.0: break 
        print printformat%(t,x,xd)

 Ein Testlauf mit obigen Angaben und einem cw-Wert von 0.43 (etwa ein Gummiball) ergibt

>>> fall(1000,0,0.43)

t = 0.00, x = 1000.0000, v =   0.0000
t = 1.00, x = 995.1530,  v =  -9.5801
t = 2.00, x = 981.2590,  v = -17.9237
t = 3.00, x = 959.9493,  v = -24.3597
t = 4.00, x = 933.1840,  v = -28.8748
t = 5.00, x = 902.7185,  v = -31.8354
t = 6.00, x = 869.8803,  v = -33.6918
t = 7.00, x = 835.5753,  v = -34.8233
t = 8.00, x = 800.3840,  v = -35.5010
t = 9.00, x = 764.6645,  v = -35.9028
t = 10.00, x = 728.6328, v = -36.1396
t = 11.00, x = 692.4177, v = -36.2785
t = 12.00, x = 656.0948, v = -36.3599
t = 13.00, x = 619.7090, v = -36.4075
t = 14.00, x = 583.2864, v = -36.4353
t = 15.00, x = 546.8422, v = -36.4516
t = 16.00, x = 510.3854, v = -36.4611
t = 17.00, x = 473.9213, v = -36.4667
t = 18.00, x = 437.4528, v = -36.4699
t = 19.00, x = 400.9819, v = -36.4718
t = 20.00, x = 364.5095, v = -36.4729
t = 21.00, x = 328.0363, v = -36.4735
t = 22.00, x = 291.5625, v = -36.4739
t = 23.00, x = 255.0885, v = -36.4741
t = 24.00, x = 218.6143, v = -36.4743
t = 25.00, x = 182.1400, v = -36.4743
t = 26.00, x = 145.6656, v = -36.4744
t = 27.00, x = 109.1912, v = -36.4744
t = 28.00, x =  72.7168, v = -36.4744
t = 29.00, x =  36.2424, v = -36.4744
t = 30.00, x =  -0.2320, v = -36.4744

Wir erkennen die verlängerte Fallzeit und das Erreichen einer konstanten Endgeschwindigkeit, wenn die Reibungskraft schließlich die Schwerkraft ausgleicht. Das ist es, was ein Fallschirmspringer ausnützt. Vielleicht kannst Du den ungefähren cw-Wert eines Fallschirms herausfinden? Du musst im Programm aber auch noch die Größe des Schirms und die Masse des Springers passend statt der Werte für unsere Kugel eintragen.

3. Schiefer Wurf

Spannend ist der schiefe Wurf. Wie wirkt sich hier die Luftreibung aus? Mit einem cw-Wert von 0 kannst Du sie zu Vergleichszwecken einfach ausschalten. Was ist der optimal Abschusswinkel für möglichst große Wurfweite - ohne Reibung sind es 45°. Ist es mit Reibung mehr oder weniger? Kugelstoßer werfen ihre Kugel nicht aus Höhe Null ab - was ist hier der optimale Winkel ohne und mit Luftreibung?

Das Programm rechnet nun zweidimensional. Das Runge-Kutta-Verfahren bekommt deshalb auch doppelt so viele Parameter für Ort und Geschwindigkeit und die Funktion f liefert zwei Komponenten zurück: die Ableitung in x- und in y-Richtung.

## schiefer Wurf mit Luftreibung

from math import *


###########################################################
# Solution of y''(t) = f(t,y,y'), y(t0)=y0, y'(t0)=y'0
# RKN4 y''=f(t,x,y,x',y')
def rkn4_1(t,x,y,xd,yd,f,h):
    k1x,k1y = f(t,x,y,xd,yd)
    k2x,k2y = f(t+h/2.0,
                x+h/2.0*xd+h*h/8.0*k1x,y+h/2.0*yd+h*h/8.0*k1y,
                xd+h/2.0*k1x,yd+h/2.0*k1y)
    k3x,k3y = f(t+h/2.0,
                x+h/2.0*xd+h*h/8.0*k1x,y+h/2.0*yd+h*h/8.0*k1y,
                xd+h/2.0*k2x,yd+h/2.0*k2y)
    k4x,k4y = f(t+h,
                x+h*xd+h*h/2.0*k3x,y+h*yd+h*h/2.0*k3y,
                xd+h*k3x,yd+h*k3y)
    xnew = x+h*xd+h*h/6.0*(k1x+k2x+k3x)
    ynew = y+h*yd+h*h/6.0*(k1y+k2y+k3y)
    xdnew = xd+h/6.0*(k1x+2.0*k2x+2.0*k3x+k4x)
    ydnew = yd+h/6.0*(k1y+2.0*k2y+2.0*k3y+k4y)
    return (t+h,xnew,ynew,xdnew,ydnew)


###########################################################
# Fall einer Kugel mit Newton'scher Reibungskraft
# F = ma , d.h. xd' = -g + cw.rho*A*v**2/m

def f(t,x,y,vx,vy):
    v = (vx*vx+vy*vy)**0.5
    ablx = -factor*v*vx    # entspricht v**2*(vx/v) 
    ably = -g -factor*v*vy 
    return ablx,ably
        


def wurf(hoehe=2,geschwindigkeit=20,winkel=55,cw=0.43):

    cw = cw         # Widerstandsbeiwert
    m = 0.8         # Masse
    r = 0.03        # Radius
    A = 4.0*pi*r**2 # Querschnittsfläche
    rho = 1.213     # Luftdichte
    global g
    g = 9.81        # Fallbeschleunigung
    global factor
    factor = cw*rho*A/m
    vx = geschwindigkeit*cos(winkel*pi/180.0)
    vy = geschwindigkeit*sin(winkel*pi/180.0)

    t = 0.0                     # laufende Zeit
    x = 0.0                     # Startwert Ort x
    y = float(hoehe)            # y
    xd = vx                     # Startwert Geschwindigkeit xd
    yd = vy                     # yd
    h = 0.01                    # Integrations-Zeitschrittweite


    printformat = "t = %6.3f, x = %7.3f, y = %7.3f, vx = %7.3f,"+
                  "vy = %7.3f, v = %7.3f"
    print printformat%(t,x,y,xd,yd,(xd**2+yd**2)**0.5)
    
    while y>=0.0:
        for i in range(100):
            t0,x0,y0,xd0,yd0 = t,x,y,xd,yd  # alte Werte merken
            t,x,y,xd,yd = rkn4_1(t,x,y,xd,yd, f,h)
            if y<=0.0: break 
        print printformat%(t,x,y,xd,yd,(xd**2+yd**2)**0.5)

    ## Aufschlag genauer bestimmen, lineare Interpolation genügt
    nullwert = y0/(y0-y)
    t1 = (t-t0)*nullwert
    # letzter Schritt nochmals mit t1 statt h
    t,x,y,xd,yd = rkn4_1(t0,x0,y0,xd0,yd0, f,t1) 
    print "\nAufschlag"
    print printformat%(t,x,y,xd,yd,(xd**2+yd**2)**0.5)

Zusätzlich habe ich hier einen abschließenden Rechenschritt angefügt, der die Aufschlagstelle unseres Balls (Höhe y=0) genauer bestimmt - durch eine lineare Interpolation ('nullwert' gibt an, bei wieviel Prozent der letzten y-Änderung der Nulldurchgang gewesen wäre, diesen Prozentsatz nehmen wir als letzte Zeitschrittweite)

Ein Beispiel: Ein Ball wird aus 2 Metern Höhe mit 20 m/s Geschwindigkeit unter 55° (gegen die Waagrechte) und obigem cw-Wert geworfen.

>>> wurf(2,20,55,0.43)
t =  0.000, x =   0.000, y =   2.000, vx =  11.472, vy =  16.383, v =  20.000
t =  1.000, x =  10.791, y =  12.684, vx =  10.246, vy =   5.308, v =  11.539
t =  2.000, x =  20.652, y =  13.008, vx =   9.499, vy =  -4.533, v =  10.525
t =  3.000, x =  29.741, y =   3.919, vx =   8.630, vy = -13.440, v =  15.972
t =  3.270, x =  32.033, y =  -0.005, vx =   8.346, vy = -15.601, v =  17.693

Aufschlag
t =  3.270, x =  32.031, y =   0.000, vx =   8.346, vy = -15.599, v =  17.691

Er ist also 3.27 Sekunden in der Luft, fliegt 32 Meter weit und kommt mit einer Geschwindigkeit von 17.69 m/s auf dem Boden auf. Vergleiche mit cw=0 !

Nun kannst Du unterschiedliche Situationen simulieren: Schlagball, Kugelstoßen, Tischtennis,...

Erweiterung:

Verpasse dem Programm eine GUI, in der man alle Parameter bequem einstellen kann und auf Knopfdruck die berechnete Bahn in einem Canvas aufgezeichnet bekommt. Eventuell sogar mehrere Bahnen übereinander, zum Vergleichen der Wurfweiten.