Das deutsche Python-Forum

Code: Alles auswählen

#gaussnormal F(x), f(x)
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
x = np.linspace(-3, 3,100)
###########
plt.plot(x, special.erf(x)) , plt.plot(x, special.erf(1/x)),plt.plot(x, special.erf(np.sin(x**2)))
##############
plt.xlabel('$x oder 1/x  $')
plt.ylabel('$erf(x,1/x)$')
plt.title("  x  und  1/x   ") 
plt.grid(True) 
plt.show()

Code: Alles auswählen

from sympy import *
x = symbols('x')


f= 1/2*x
print("f(x) =", f)
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print(log(x) )


f = 1+1/x

print("f(x) =", f)
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print(log(x) )


f=  3/x
print("f(x) =", f)
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print(log(x) )

Code: Alles auswählen

#  "Exotische Funktionen"  werden
#   1.  integriert
#   2.  abgeleitet
#  "Exotische Funktionen"  werden
#   1.  integriert
#   2.  abgeleitet


from sympy import *
import numpy as np
x = symbols('x')

f = 1/x 
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print()               




print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print()               

f=   3*np.e **cos(x**2) -( np.pi**(- 2*x) + np.pi + sin(x))
                           
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))


from sympy import *
import numpy as np
x = symbols('x')

f = 1/x 
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print()               



print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print()               

f=   3*np.e **cos(x**2) -( np.pi**(- 2*x) + np.pi + sin(x))
                           
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))

Code: Alles auswählen

#   
#  "Exotische Funktionen"  werden
#   1.  integriert
#   2.  abgeleitet


from sympy import *
import numpy as np
x = symbols('x')

f = 1/x 
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print()               

f=  (x**2 + 3*np.sin(np.e**2)  /  (3*x  + np.e**cos(np.pi)))

print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))
print()               

f=   3*np.e **cos(x**2) -( np.pi**(- 2*x) + np.pi + sin(x))
                           
print("Stammfunktion: F(x) =", integrate(f, x))
print("Ableitung: f'(x) =", diff(f, x))

Code: Alles auswählen

# Zum Einstieg in PLotly diesen Code    ausführen
import plotly.express as px 
# Creating the Figure instance
fig = px.line(x=[1,2, 3], y=[1, 2, 3]) 
# printing the figure instance
print(fig)
#Jetzt steht  PLOTLY  offen für  Python-Anwender

#gleich ein Beispiel
import plotly.express as px 

# Data to be plotted
df = px.data.iris() 
# Plotting the figure
fig = px.scatter_3d(df, x = 'sepal_width', 
                    y = 'sepal_length', 
                    z = 'petal_width', 
                    color = 'species') 
 
fig.show()
 

Code: Alles auswählen

import plotly.express  as  px

z  =  [[ .1 ,  .3 ,  .5 ,  .7 ,  .9 ], 
     [ 1 ,  .8 ,  .6 ,  .4 ,  .2 ], 
     [ .2 ,  0 ,  .5 ,  .7 ,  .9 ], 
     [ .9 ,  .8 ,  .4 ,  .2 ,  0 ], 
     [ .3 ,  .4 ,  .5 ,  .7 ,  1 ]]

fig  =  px.imshow ( z , text_auto = True )
fig.show()

Code: Alles auswählen

mport numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
import scipy.integrate as integrate

#Schritt 2: Wir definieren die Parameter.

g = 9.81 # [m/s^2]
l = 0.5 # [m]
gamma = 0. # [1/s] damping coefficient
omega0 = np.sqrt(g/l)

#Schritt 3: Wir definieren die Anfangsbedingungen, d. h.
# Anfangswerte für Winkel thetaund Winkelgeschwindigkeit omega.

theta= 45 # [degree]  # initial angle
theta = theta*np.pi/180
omega = 0  # initial angular velocity

#Schritt 4: Wir möchten uns ein Bild von der Zeitskala als Funktion des Anfangswinkels machen.
#Dazu berechnen wir die Schwingungsdauer Tperdurch Quadratur. #
#Mit diesem Wert legen wir dann die maximale Zeit Tmaxder Simulation als Vielfaches von fest Tper.

def f(x):
    return 1/np.sqrt(np.cos(x)-np.cos(theta))

Tper = integrate.quad(f ,0,theta )
T0 = 2.0*np.sqrt(2.0)*Tper[0]/omega0
Tmax = 8*T0
print('Tmax=',Tmax)

#Schritt 5: Wir definieren grafische und Simulationsparameter.

xmax =  1.05*l
ymax =  1.05*l
Nstep= 3600   # number of time steps / number of frames
dt= Tmax/Nstep  # time increment
d=0.05  # size of fixed support

#Schritt 6: Wir legen die Figur, die Achse und das Linienobjekt fest, das wir animieren möchten. Dieses Linienobjekt ist vorerst leer. #

#Wir zeichnen außerdem ein rechteckiges Objekt der Größe 2d x 2d, das am Ursprung zentriert ist.

fig = plt.figure()
ax = plt.axes(xlim=(-xmax, xmax), ylim=(-ymax, ymax))

line, = ax.plot([], [], 'o', ls='-', ms=8, markevery=[0,-1])
xR,yR = -d,-d # lower left corner of rectangle
dx,dy= 2*d, 2*d  # dimensions of rectangle
rect = plt.Rectangle((xR,yR), dx,dy, facecolor= 'black', edgecolor= 'black')
ax.add_patch(rect)
plt.axis('off')
ax.set_aspect('equal')

#Schritt 7: Wir definieren eine Funktion, die die Bewegungsgleichungen modelliert

def F(X,t):
    [theta, omega] = X
    dthetadt = omega
    domegadt = -g/l* np.sin(theta) - gamma *omega
    return dthetadt, domegadt

#Schritt 8: Wir definieren eine Funktion, um die Initialisierung des Plots jedes Frames einzurichten.

def init():
    line.set_data([], [])
    return line,

#Schritt 9: Wir definieren die Animationsfunktion. Bei jedem Wert der Ganzzahl iberechnen wir die Pendelposition durch Integration der ODE von Zeit timeIn = i*dtzu Zeit timeOut= (i+1)*dt. Dazu rufen wir scipyden ODE-Integrator auf odeint. Sobald wir neue Werte der Variablen thetaund haben omega, können wir die neuen kartesischen Koordinaten des Pendels berechnen und das Linienobjekt aktualisieren.

def animate(i):
    
    global theta, omega
    timeIn= i*dt
    timeOut= (i+1)*dt
    Xinit = [theta, omega]
    odeout = integrate.odeint(F, Xinit, [timeIn, timeOut])[1]
    theta = odeout[0] 
    omega = odeout[1]
    theta = (theta + np.pi) % (2*np.pi) - np.pi # bring theta back to [-pi, pi]
    x,y = l*np.sin(theta), -l*np.cos(theta)
    line.set_data([0, x], [0,y]) # modify the pendulum's position    
    return line,

#Schritt 10: Wir können nun den Animator mit der blitauf gesetzten Variable aufrufen True. Dadurch wird das animObjekt berechnet. Anschließend speichern wir das Video in einer MP4-Datei. Dieses Video können wir in unserem Browser ansehen.

anim = animation.FuncAnimation(fig, animate, init_func=init,frames=Nstep, interval=20, blit=True)
plt.show()

Code: Alles auswählen

Projektilbewegung

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
#Schritt 2: Wir definieren die notwendigen Parameter, wie
#Wir erhöhen die Zeit mit einer Schrittweite dtentsprechend Nstep.
#Nstepdefiniert auch die Anzahl der Frames der Animation.

v0 = 50 # [m/s]
g = 9.81 # [m/s^2]
alpha= 50 # [deg]
D= 220 # [m]
d= 10 # [m]
alpha = alpha*np.pi/180 # Winkel
Tmax = D/(v0*np.cos(alpha)) # Dauer Flugbahn
xmax =  D+2*d #  2d 
ymax =  0.55*(v0*np.sin(alpha))**2 /g                  
Nstep= 400   #  Anzahl Zeitschritte /number of frames
dt= Tmax/Nstep  
h= D*np.tan(alpha)-0.5*g*(D/v0/np.cos(alpha))**2 # height of rectangular object
#Schritt 3: Wir richten die Figur, die Achse und das zu animierende Plot-Element ein.
#Dieses Element ist ein leeres Linienobjekt mit den gewünschten Attributen
#Das Linienobjekt wird während der Simulation zu einem späteren Zeitpunkt aktualisiert (Schritt 6).
#Wir fügen dem Plot auch das statische rechteckige Objekt hinzu.

fig = plt.figure()
ax = plt.axes(xlim=(0, xmax), ylim=(0, ymax))
plt.xlabel('x(m): horizontal distance')
plt.ylabel('y(m): height')
line, = ax.plot([], [], 'o', ls='--', ms=8, markevery=[0,-1])
xR,yR = D-d,0. # lower left corner of rectangle
dx,dy= 2*d, h # dimensions of rectangle
rect = plt.Rectangle((xR,yR), dx,dy, facecolor= 'red', edgecolor= 'red')
ax.add_patch(rect)
#Schritt 4: Wir definieren eine Initialisierungsfunktion, um den Hintergrund jedes Frames zu plotten.
#Diese Funktion setzt das Linienobjekt mit leeren Daten. #
#Durch die Rückgabe des Linienobjekts weiß der Animator nun,
#welches Plotobjekt während der Animation aktualisiert werden muss.

def init():
    line.set_data([], [])
    return line,

#Schritt 5: Wir definieren die Animationsfunktion, die das Linienobjekt mit der Zeit aktualisiert.
#Sie hat ein einziges Argument: die Ganzzahl i, die die Frame-Nummer definiert.
#Die Zeit wird anhand dieser Ganzzahl definiert i, gemäß einer bestimmten Skala. 
#bIn unserem Fall wird das Linienobjekt durch 2 Arrays definiert
#und parametrisiert mit der Zeit gemäß den Gleichungen (1-2)
# Diese Funktion wird vom Animator sequenziell aufgerufen.

def animate(i):
    t  = np.linspace(0, i*dt, 2*i)
    x = v0*np.cos(alpha)*t
    y = v0*np.sin(alpha)*t -0.5*g*t**2
    line.set_data(x, y)
    return line,
#Schritt 6: Jetzt rufen wir den Animator auf FuncAnimation.
#Das resultierende Objekt muss einer Variablen zugewiesen werden, hier anim.
#Wir definieren die Animation durch die Anzahl der Frames Nstep  #mit einer Verzögerung von 20 ms zwischen den Frames. Das blitSchlüsselwort ist wichtig: Es weist den Animator an, nur die Teile der Handlung neu zu zeichnen, die sich geändert haben. Durch die Einstellung blit=Truewerden die Animationen viel schneller angezeigt und die resultierende Videodatei wird kleiner. Der optionale saveBefehl speichert die Daten in einer MP4-Videodatei. Der showBefehl zeigt das Endergebnis an.

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=Nstep, interval=20, blit=True)
plt.grid()
plt.title("Bahn des Projektils zum Zielbereich")
plt.show()

Code: Alles auswählen

print("Das aktuelle Datum ")
import time
print(time.strftime("%d.%m.%Y %H:%M:%S"))                     #####  tt,mm,jj,sec
now = time.localtime()                                                             #####  tt,m,JJ
print( now.tm_mday,".",now.tm_mon,".", now.tm_year   )


import datetime
from datetime import  timedelta
now = datetime.datetime.now()
now
print("Ausgangsdatum, also heute: " ,now)
print()
 
print("Nach 7 Tagen:",now + timedelta(days=7))
print(now + timedelta(days=7))
print()
print("Vor  12 Tagen:",now - timedelta(days=12)) 
print(now  - timedelta(days=12))
print()
print("Vor  712 Tagen:",now - timedelta(days=12)) 
print(now -  timedelta(days= 712))

Code: Alles auswählen

import datetime
Zeit = datetime.datetime.now()
print(Zeit)
Datum = Zeit.strftime('%d.%m.%Y')             #  US Format        
print(Datum)
Uhrzeit = Zeit.strftime('%H:%M:%S')           #    Dt.Format  
print(Uhrzeit)
Genau = Zeit.strftime('%d:%m.%Y:%H:%M:%S')            #Dt.Format
print(Genau)
 

Code: Alles auswählen

print("Das heutige Datum    - zunächst  US-Format")
from datetime import date

today = date.today()
print(today)
print()
print("Nach Konvertierung  in EU/Dt - Format ")
print("                            1.   Tag  Monat  Jahr")

d1 = today.strftime("%d/%m/%Y")
print("d1 zeigt     =", d1)
print()
print("                             2.   Monat Tag  Jahr")
# Textual month, day and year	

d2 = today.strftime("%B %d, %Y")
print("d2 =", d2)
print()

print("                             3.   Monat  Tag   Jahr")
d3 = today.strftime("%m/%d/%y")
print("d3 =", d3)
print()

print("(abgekürzt)           4.   Monat   Tag   Jahr")
# Month abbreviation, day and year	
d4 = today.strftime("%b-%d-%Y")
print("d4 =", d4)

Code: Alles auswählen

Folgende Tabelle enthält sind alle Formatcodes aufgeführt, die man verwenden kann.

Richtlinie	Bedeutung	Beispiel
%a	Abgekürzter Wochentagsname.	So, Mo, ...
%A	Vollständiger Wochentagsname.	Sonntag, Montag, ...
%w	Wochentag als Dezimalzahl.	0, 1, ..., 6
%d	Tag des Monats als mit Nullen aufgefüllte Dezimalzahl.	01, 02, ..., 31
%-d	Tag des Monats als Dezimalzahl.	1, 2, ..., 30
%b	Abgekürzter Monatsname.	Jan., Feb., ..., Dez.
%B	Vollständiger Monatsname.	Januar, Februar, ...
%m	Monat als mit Nullen aufgefüllte Dezimalzahl.	01, 02, ..., 12
%-m	Monat als Dezimalzahl.	1, 2, ..., 12
%y	Jahr ohne Jahrhundert als mit Nullen aufgefüllte Dezimalzahl.	00, 01, ..., 99
%-y	Jahr ohne Jahrhundert als Dezimalzahl.	0, 1, ..., 99
%Y	Jahr mit Jahrhundert als Dezimalzahl.	2013, 2019 usw.
%H	Stunde (24-Stunden-Uhr) als mit Nullen aufgefüllte Dezimalzahl.	00, 01, ..., 23
%-H	Stunde (24-Stunden-Uhr) als Dezimalzahl.	0, 1, ..., 23
%I	Stunde (12-Stunden-Uhr) als mit Nullen aufgefüllte Dezimalzahl.	01, 02, ..., 12
%-I	Stunde (12-Stunden-Uhr) als Dezimalzahl.	1, 2, ... 12
%p	Ortsangabe: AM oder PM.	Vormittags, nachmittags
%M	Minute als mit Nullen aufgefüllte Dezimalzahl.	00, 01, ..., 59
%-M	Minute als Dezimalzahl.	0, 1, ..., 59
%S	Sekunde als mit Nullen aufgefüllte Dezimalzahl.	00, 01, ..., 59
%-S	Sekunde als Dezimalzahl.	0, 1, ..., 59
%f	Mikrosekunde als Dezimalzahl, links mit Nullen aufgefüllt.	000000 - 999999
%z	UTC-Offset im Format +HHMM oder -HHMM.	 
%Z	Name der Zeitzone.	 
%j	Tag des Jahres als mit Nullen aufgefüllte Dezimalzahl.	001, 002, ..., 366
%-j	Tag des Jahres als Dezimalzahl.	1, 2, ..., 366
%U	Wochennummer des Jahres (Sonntag als erster Tag der Woche). Alle Tage im neuen Jahr vor dem ersten Sonntag werden als Woche 0 betrachtet.	00, 01, ..., 53
%W	Wochennummer des Jahres (Montag als erster Tag der Woche). Alle Tage im neuen Jahr vor dem ersten Montag werden als Woche 0 betrachtet.	00, 01, ..., 53
%c	Entsprechende Datums- und Zeitdarstellung des Gebietsschemas.	Montag, 30. September 2013, 07:06:05 Uhr
%x	Entsprechende Datumsdarstellung des Gebietsschemas.	30.09.13
%X	Entsprechende Zeitdarstellung des Gebietsschemas.	07:06:05
%%	Ein wörtliches „%“-Zeichen.	%   usw.

Code: Alles auswählen

import datetime

current_date = datetime.datetime.now()
current_day = current_date.strftime("%A")
print("Today is", current_day)

current_date = datetime.datetime.now()
current_day = current_date.strftime("%Y")
print("Today is", current_day)

current_date = datetime.datetime.now()
current_day = current_date.strftime("%d.%m.%Y")
print("Today is", current_day)

current_date = datetime.datetime.now()
current_day = current_date.strftime("%H.%Y")
print("Today is", current_day)

Code: Alles auswählen

A.#Programm  sucht und findet Wörter und zählt diese
#Die ausgewählte  - datei.py  -  muss im gleichen Ordner oder per Pfad erreichbar sein.
#                              z.B. den oben gezeigten Text      Schiller.py   eingeben
def suchen_finden():
    """words_suchen_finden() - zählt alle Wörter."""

    file_name = input("Namen der Datei eingeben: ")
    # Datei öffnen  und in Liste speichern
    with open(file_name, "r") as f:
        word_list = f.read().split()
    # Anzeige des Inhaltes der Python-Datei 
    suchen_finden_dict = {}
    for word in word_list:
        # Vorkommen auszählen
        suchen_finden_dict[word] = suchen_finden_dict.get(word, 0) + 1
    # Ergebnisse
    print(
        f"File {file_name} has {len(word_list)} words, "
        f"{len(suchen_finden_dict)} are unique:"
    )
    print(suchen_finden_dict)

if __name__ == "__main__":
       suchen_finden() 

Code: Alles auswählen

Programm   zählt  alle  Wörter 
Friedrich Schiller: Der Handschuh 
Vor seinem Löwengarten
Das Kampfspiel zu erwarten
Saß König Franz
Und um ihn die Großen der Krone  
Und rings auf hohem Balkone 
Die Damen in schönem Kranz   
Und wie er winkt mit dem Finger
Auftut sich der weite Zwinger
Und hinein mit bedächtigem Schritt
Ein Löwe tritt
Und sieht sich stumm
Rings um
Mit langem Gähnen,
Und schüttelt die Mähnen
Und streckt die Glieder
Und legt sich nieder.

Code: Alles auswählen

                                     Simulation   einer   Integration
#                                        am     Beispiel   von  sin(x)
#                                      Ergebnis             Genauigkeit
import scipy                    
import random 
import numpy as np 
  
# Grenze der Integration 
a = 0
b = np.pi                        # Grenzen 0 bis  = Pi
#Einzelschritte
N = 1000                       #  mit unterschiedlichem   N   Ergebnis vergleichen
  
# array of zeros of length N 
ar = np.zeros(N)   
 # Iteration  über N         #Grenzen   a un b                                           
for i in range (len(ar)): 
    ar[i] = random.uniform(a,b) 
  
# variable to store sum of the functions of  
# different values of x 
integral = 0.0
  
# function to calculate the sin of a particular 
# value of x 
def f(x): 
    return np.sin(x) 
  
# iterates and sums up values of different functions 
# of x 
for i in ar: 
    integral += f(i) 
  
# we get the answer by the formula derived adobe 
ans = (b-a)/float(N)*integral 
  
# prints the solution 
print ("Wert mit  Simulation   berechnet    :      {}."   .format(ans)) 
print ("mit           Integration : ", 2.000000000)
print("Genauigkeit steigt mit    N")

Code: Alles auswählen

#                                 Monte-Carlo-Simulationen  (aus   der Stochastik )

#                                  Das Roulette Rad hat
#                                              18 rote,
#                                    18 schwarze Fächer
#                                              1   grünes Fach,
#                                   also insgesamt 37 Fächer

#  Wenn wir    beim Roulette immer          Black     wählen,      wie oft würden wir gewinnen?“


import random

def play_roulette():
    total_slots = 37
    red_probability   = (18 / total_slots) * 100
    black_probability = (18 / total_slots) * 100
    green_probability = (1 / total_slots) * 100
    
    possible_outcomes = ["red", "black", "green"]
    probabilities = [red_probability, black_probability, green_probability]
    
    outcome = random.choices(
        possible_outcomes,
        weights=probabilities,
        k=1
    )[0]
    
    return outcome
    
#HIer begInnt die Simulation    
def perform_simulation(n_times=1, choice="black"):   
    results = { "red": 0, "black": 0, "green": 0 }
    
    for i in range(n_times):
        outcome = play_roulette()
        results[outcome] += 1
           
    win_percentage = results[choice] / n_times
        
    return results, win_percentage

        
if __name__ == "__main__":
    results, win_percentage = perform_simulation(n_times=1000,   choice="black")
                                                      
    print(results)
    print(" In Prozent  :",win_percentage)
 

Code: Alles auswählen

#                SIMULATION    mit    DGL   als  'Tool'
import numpy as np
from scipy.integrate import odeint
#Patienten,Infektionen, Geheilte

import matplotlib.pyplot as plt
#Die Anfangsbedingungen werden  zuerst gesetzt: 
S0, I0, R0 = 500,1, 1
#und ebenso die Modellparameter alpha und beta:

alpha, beta  = 0.001, 0.01
#Die rechten Seiten der Differentialgleichungen
# werden als eine vektor-wertige Funktion geschrieben:

def deriv(u, t, alpha, beta):
   S, I, R = u
   dSdt = -alpha * S * I
   dIdt = alpha * S * I - beta  * I
   dRdt = beta * I
   return dSdt, dIdt, dRdt
#dann setzen wir auch die Anfangsbedingungen als Vektor der Werte S0, I0 und R0:

y0 = S0, I0, R0
#wir möchten die Lösung zu diskreten Zeitschritten bestimmen und 
#legen dafür ein Zeitintervall und Zeitschritte fest:

T = np.linspace(0, 100, 200)
#Anschließend lösen wir das S-I-R Modells zu den diskreten Zeitpunkten:
ret = odeint(deriv, y0, T, args=(alpha, beta))
S, I, R = ret.T
# Jetzt wird die Anzahl für jeden Zeitpunkt bestimmt 
 

##Mit Hilfe der matplotlib Erweiterung können wir das Ergebnis visualisieren:

plots = plt.plot(T, S, 'b', alpha=0.5, lw=2, label='Probanden') \
   + plt.plot(T, I, 'r', alpha=0.5, lw=2, label='infiziert') \
   + plt.plot(T, R, 'g', alpha=0.5, lw=2, label='geheilt,immun')
plt.xlabel('Zeit [Tage]')
plt.ylabel('Probanten')
plt.legend(handles=plots)
plt.grid()
plt.title("Ist  DGL als 'Tool  für Simulation'  geeignet ?   JA") 
plt.show()

Code: Alles auswählen

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
 
# das Differentialgleichungssystem mit dem Vektor U=(u1, u2)
# die Funktion gibt einen Vektor mit den Ableitungen zurück.
def DGLSystem(U, y):
    u1, u2 = U
    mu = 5.0                  
    dpdx = -200.0
    du1dy = u2
    du2dy = dpdx / mu
    return [du1dy, du2dy]

# Randwerte (gegeben)
u1_0 = 0.0
u1_H = 0.0
H    = 0.1

# Randwert (geschätzt)
u2_0 = 2.0

# diskrete y-Werte (hier 50 über die Spalthöhe)
yi = np.linspace(0, H, num=50)

# Lösen des DGL-Systems mit  scipy odeint
 
U = odeint(DGLSystem, [u1_0, u2_0], yi)

# Verlauf von u_1 darstellen
 
plt.figure(figsize=(6, 8))
plt.plot(U[:,0], yi)
plt.plot([u1_H], [H], 'ro')
plt.xlabel('$u_1$')
plt.ylabel('$y$')
plt.grid()
plt.title("Roter Punkt == Zielpunkt")
plt.show();

#Die Schätzung des Randwerts $u_2(0)$ lässt sich natürlich auch automatisieren.
#Eine sehr einfache Möglichkeit ist die Anpassung des Randwertes in Abhängigkeit 
#von der berechneten Abweichung:

 
# Vorgegebene Toleranz
tol = 0.0001
faktor = 1.0
iterationen = 0

# Randwerte (gegeben)
u1_0 = 0.0
u1_H = 0.0
H    = 0.1

# Randwert  probieren: 1 , 2 oder  3 
u2_0 = 3

# Lösen des gewöhnlichen DGL-Systems 
#  mit (SciPy odeint
 
fehler = 2 * tol

plt.figure(figsize=(6, 8))

while abs(fehler) > tol and iterationen < 100:
    U = odeint(DGLSystem, [u1_0, u2_0], yi)
    plt.plot(U[:,0], yi)

    fehler = (u1_H - U[-1,0])
    print(fehler)
    u2_0 = u2_0 + fehler * faktor
    iterationen += 1

print ("Bis zum Ziel wurden ", iterationen, " Iterationen benötigt.")
plt.title("Erste Annäherung an das Ziel")
# Verlauf von u_1 darstellen
plt.plot(U[:,0], yi)
plt.plot([u1_H], [H], 'ro')
plt.xlabel('$u_1$')
plt.ylabel('$y$')
plt.grid() 
plt.title("Der rote Punkt ist das Ziel")

plt.show()

Code: Alles auswählen

# Simulation Freier Fall

g =  9.81 # m/s² (Fallbeschleunigung Erde)
## Anfangsbedingungen:
s_0 =  200    # m    (Anfangshöhe)
print("Fallhöhe in  Meter :",s_0)
v_0 =      0    # m/s  (Anfangsgeschwindigkeit)
## Genauigkeit:
dt =       0.1  # s    (Zeitschritt)
## Simulation:
a = -g
s =  s_0
v =  v_0
t =  0
while s > 0:
    s += v * dt                    # neues s aus dem alten Wert von v berechnen
    v += a * dt                    # neues v aus dem alten Wert von a berechnen
    t += dt
## Textausgabe:
print("Falldauer: ", round(t,1), " s")
print("Endgeschwindigkeit: ", round(v,), "m/s")

 ## Simulation:
list_t = []            
list_s = []
a = -g
s = s_0
v = v_0
t = 0
while s > 0:
    list_t.append(t)
    list_s.append(s)
    s += v * dt                   
    v += a * dt                    
    t += dt
  

## Diagramm:
import matplotlib.pyplot as plt
plt.plot(list_t, list_s, 'b.')
plt.ylabel('Höhe (m)')
plt.xlabel('Zeit (s)')
plt.title('Simulation: Freier Fall aus '+str(s_0)+'m')
plt.grid()
plt.show()