kleine Mathe-Spielereien

[OSWALD]

4.03.2023
Ich habe es gewagt, mich in die Thematik
'Differemtialgleichungen '
zu wagen und bin nun soweit gekommen,
dass ich glaube das Prinzip 'Lösung von Dgl mit Hilfe von scipy u. numpy'
in etwa vertanden zu haben. Aus scipy habe ich den ersten Code *Pendelbewegung'
entnommen.
Es ist eine einfache Dgl, die auch sehr gut erklärt ist.
Ich habe weitere einfache Dgl 's studiert und werde auch diese
als Beispiele vorstellen . Mit dem ''Werkzeug' Differential- und
Integralrechnung' sollten die math. technischen Voraussetzungen
für eine erffolgreiche Einarbeitung gegeben sein.
Gute Zeit OSWALD

Code: Alles auswählen

#theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0
#wobei b und c positive Konstanten sind 
#und ein Strich (') eine Ableitung bezeichnet. 
#Um diese Gleichung mit  odeint zu lösen 
#müssen wir sie zunächst in ein System von Gleichungen erster Ordnung umwandeln. 
#Durch Definition der Winkelgeschwindigkeit erhalten wir
#                  das System:       omega(t) = theta'(t)
#                                    theta'(t) = omega(t)

#                                    omega'(t) = -b*omega(t) - c*sin(theta(t))
#                  der Vektor sei   [ theta , omega ].  

import numpy as np
def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

#Wir nehmen an, dass die Konstanten b = 0,25 und c = 5,0 sind:

b = 0.25
c = 5.0
#Für die Anfangsbedingungen nehmen wir an, 
#dass das Pendel mit Theta(0) = pi – 0,1 nahezu vertikal steht 
#und anfänglich in Ruhe ist, 
#     also omega(0) = 0. Dann ist der Vektor der Anfangsbedingungen
y0 = [np.pi - 1.0,0]

#Wir werden eine Lösung mit 101 gleichmäßig beabstandeten Abtastungen 
#im Intervall 0 <= t <= 10 generieren. Unsere Zeitreihe ist also:

t = np.linspace(0, 10, 101)

#Rufen Sie   ++***odeint ***   auf , um die Lösung zu generieren. 
#Um die Parameter      b und c an pend       zu übergeben , 
#geben wir sie mit      odeintdem      Argument args an .

from scipy.integrate import odeint
sol = odeint(pend, y0, t, args=(b, c))
#Die Lösung ist ein Array mit der Form (101, 2).
# Die erste Spalte ist Theta(t) und die zweite ist Omega(t) . 
# Der folgende Code zeichnet beide Komponenten.

import matplotlib.pyplot as plt
plt.plot(t, sol[:, 0], 'b', label='theta(t)')
plt.plot(t, sol[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel(" omega(t) = theta'(t)  ") 
plt.grid()
plt.show()

[OSWALD]

6.03.2023

Hier zeige ich eine der einfachsten Dgl,s , die ich kenne.
y = meine_Funktion , y0, t)
y soll der Weg sein, hier ymax = 5 in m
k = Anfangwert. Wenn k =0, dann bei Y0 sehen wir nur eine Gerade
je größer k, umso steiler die Kurve
t = hier die Zeit
' meine_Fiuktion ' ist also eine Wachstumsfunktion
sie hat natürlich auch eine Umkehrfunktion
dt = der Zeitschritt
die function t = dy/dt
Gute Zeit OSWALD

Code: Alles auswählen

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def meine_Funktion(y,t):
    dt= 1
    k = 2.3
    dydt = k * y
    return dydt

# initial condition
y0 = 5

# time points
t = np.linspace(0,2) 

# solve ODE
y = odeint(meine_Funktion,y0,t)
print(y)
# plot results
plt.plot(t,y)
plt.xlabel('Zeit')
plt.ylabel('Weg(t)')
plt.grid()
plt.show()

[OSWALD]

6.03.2023
Analytische Lösung der Dgl y'' = -g . g = 9.81 m/sec²
g = 9.81
Freier Fall im Vakuum. . Unabhängig von der Masse.
Die Geschwindigkeit ist die erste Ableitung des Ortes nach der Zeit.
Die Beschleunigung ist die erste Zeitableitung der Geschwindigkeit bzw.
die zweite Zeitableitung des Wegs. (von 0 - --->
Die Beschleunigung ist die erste Zeitableitung der Geschwindigkeit bzw.
die zweite Zeitableitung die des Ortes.
# g = dv/dt
g= d²y/dt²
Um (t) zu berechnen ,müssen die beiden Funktionen getrennt werden:
g= dv/dt und g*y/dt
Jetzt integrieren wir über die beiden Variablen und erhalten.
∫dv = ∫g*dt
v - v =g*t
v(t) =v0 +g*t

Um auch die Ortsfunktion zu bekommen,
geht man analog wie bei der Geschwindigkeitsfunktion vor.
y = dy/dt -> dy v*dt
∫dy = ∫v*dt
y(t) = y0 +v0*t +1/2*g*t*²

Daraus ergeben sich letzten Endes , die drei unten gezeigten Funktionen,
(Probleme mit der Übertragung)
############################################################
Gute Zeit OSWALD

Code: Alles auswählen

[/
import numpy as np
import matplotlib.pyplot as plt
 
 
g = 9.81                     # Erdbeschleunigung (m/s^2)
dt = 0.1                     # Zeitschritt (y)
# Anfangsbedingungen:
t = 0          
y = 0            # Fallweg (m)
v = 0            # Geschwindikeit (m/sec)
while t< 10 :
  
 v = v - g* dt
 y = y + v* dt
 t = t + dt

print("Fallhöhe ="  ,round(y ),"m")
print("Fallzeit = :", round(t),"sec")
print("Endgeschw. = :",round(v*3.6,1),"km/h")
 
############################################ 












code]

[OSWALD]

8.003.2023
Drei von mir vorgegebene und sehr einfache
Dgl-Gleichungen
stelle ich hier vor, numerisch gelöst..
Ich bin mir aber nicht unbedingt sicher,
dass die Lösungen auch korrekt sind.
Gute Zeit Oswald

Code: Alles auswählen

#Numerische Lösung von Dgl-Gleichungen

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from  numpy import *
def returns_dydt(y,t):
    dydt =  -y * t + np.cos(y)
    #dydt = -y * t + np.exp(y)
    #dydt =  y * t + np.sin(2*y)
    return dydt
  
 
y0 = 1      #Anfangswert
  
# Zeitwerte
t = np.linspace(0,15)
  
# Lösung der Dgl
y = odeint(returns_dydt, y0, t)
 
# plotten der Resultate: dt, y0, t)
print(y)  
plt.plot(t,y)
plt.xlabel("Zeit")
plt.ylabel("Y")
plt.title(" 3 verschuedene Dgl- Funktionen  zur Wahl"  )
plt.grid()
plt.show()

[OSWALD]

9.03.2023
#Hier meine Beschreibung der vorliegenden Dgl, so wie
#ich das Thema verstanden habe. Falls diese nicht richtig ist
#bitte Korrektur.

#Bei der vorliegenden Dgl handelt es sich im Grunde um
#einen Vergleich zwischen einer numerischen und einer
#analytischen Lösung.
#Für analytische Lösngen gibt es zahlreiche Verfahren, wie z.B. nach Euler oder
#das Runge Kutta-Verahren, wobei das letztere genauer sein soll.
#Konkret hier wird eine Lösung zwischen
#Scipy und RFunge Kutta verglichen.

#Vor dem Vergleich muss te die _Normal_Dgl in die analytische Form nach
#nach Runge Kutta gebracht werdn.

#Info:
# Beide Anfangswerte auf 0 setzen.
Beide derzeit ungleich, um beide Funktionen besser erkenn zu können,

Code: Alles auswählen

##########################################
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp


# Normal-Dgl gelöst mit 
ode_fn = lambda t, x: np.sin(t) + 3. * np.cos(2. * t) - x


#analytische Lösung der Gleichung
an_sol = lambda t : (1./2.) * np.sin(t) - (1./2.) * np.cos(t) + \
                    (3./5.) * np.cos(2.*t) + (6./5.) * np.sin(2.*t) - \
                    (1./10.) * np.exp(-t)
t_begin= 1  #  Anfangswert    
t_end=10.  #
t_nsamples=100
t_space = np.linspace(t_begin, t_end, t_nsamples)

x_init = 0   #Anfangswert  x. wenn beide AW gleich völlige Deckung der Funktionen

x_an_sol = an_sol(t_space) 
 
print(x_an_sol)
 
 
method = 'RK45' #available methods: 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'
num_sol = solve_ivp(ode_fn, [t_begin, t_end], [x_init], method=method, dense_output=True)

print(num_sol)


x_num_sol = num_sol.sol(t_space).T

plt.figure()
plt.plot(t_space, x_an_sol, '--', linewidth=3, label='analytisch')
plt.plot(t_space, x_num_sol, linewidth=1, label='numerisch')
plt.title('Dgl 1.Ordnung mit IVP gelöst und SciPy mit RK45 sind  exakt gleich')
plt.xlabel('t')
plt.ylabel('x')
plt.legend()
plt.grid()
plt.show()

[OSWALD]

12.3.2023
SIR-Modelle (susceptible-infected-removed model) sind beim
Einsatz bzw. Verwendung von DGL von ausserordentlicher Bedeutung.
Mit diesem selbserklaerenden SIR-Beispiel konnte ich mein Verstaendnis
füer die Loesung von Dgls mit scipy u. numpy weiter vertiefen.
Ich muss aber glechzeitig versuchen , auch analytische
Loesungen wirklich zu verstehen. Nicht jeder wird ein guter Koch ,
wenn er Kochbuecher lesen kann.

Dieser link ist 'Besessene':
https://itp.uni-frankfurt.de/~hanauske/ ... ading.html

Gruss OSWALD

'Wer immer strebend sich bemüht, den koennen wir erloesen'(FaustII,G)

[OSWALD]

hoffentlich klappt es jetzt

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
tmax=200
#Anfangswerte
S0=9970     #nicht immune Gesunde
I0=3       #Infizierte
R0=0       #Genesene
N=S0+I0+R0 #Population
#Parameter
b=0.4  #Infektionsrate
g=0.04 #Genesungsrate
#DGL-System
def dgl(t,ya):
    S,I,R=ya
    dS_dt=-b*S*I/N    #nicht immune Gesunde
    dI_dt=b*S*I/N-g*I #Infizierte
    dR_dt=g*I         #Genesene
    return [dS_dt,dI_dt,dR_dt]
#Anfangswerte
y0 = [S0,I0,R0]
t = np.linspace(0, tmax, 500)
z=solve_ivp(dgl,[0,tmax],y0,dense_output=True) 
S, I, R = z.sol(t) #1
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(t, S,'b-', lw=2, label="nicht immune Erwachsene")
ax.plot(t, I,'r--',lw=2, label="Infizierte")
ax.plot(t, R,'g:', lw=2, label="Genesene")
ax.legend(loc="best")
ax.set_xlabel("Zeit/ Tage")
ax.set_ylabel("Erwachsene")
ax.grid(True)
plt.show()        

[OSWALD]

12.3.2023
Um das Ganze etwas anschaulicher zu machen,
habe ich einmal konkrete Zahlen ausgedruckt.
Dazu habe ich die Anzahl der Probanden auf 1000
reduziert , um eine klare Darstellung zu ermöglichen.
Für jeden 50.Tag, also 20 mal werden die die Zahlen
zum Zeitpunkt X angezeigt.
Einfach einmal Variable und und andere Parameter
verändern und sehen was passiert.
So geht' ja auch in der Praxiks.
Gute Zeit OSWALD

Code: Alles auswählen

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
tmax=100
#Anfangswerte
S0=1000     #nicht immune Gesunde
I0=3       #Infizierte
R0=0       #Genesene
N=S0+I0+R0 #Population
#Parameter
b=0.4  #Infektionsrate
g=0.04 #Genesungsrate
#DGL-System
def dgl(t,ya):
    S,I,R=ya
    dS_dt=-b*S*I/N    #nicht immune Gesunde
    dI_dt=b*S*I/N-g*I #Infizierte
    dR_dt=g*I         #Genesene
    return [dS_dt,dI_dt,dR_dt]
#Anfangswerte
y0 = [S0,I0,R0]

 
t = np.linspace(0, tmax, 20) #Zählung je nach 20 Tagen
z=solve_ivp(dgl,[0,tmax],y0,dense_output=True) 
S, I, R = z.sol(t) #1
# Gezählt wird jeweils nach 2q0 Tagen  
print("Anzahl der Erwachsenen")
print(np.floor(S))
print() 
print("Anzahl der Infizierten")
print(np.floor(I))
print()
print("Anzahl der Genesenen")
print(np.floor(R))
 


fig, ax = plt.subplots(figsize=(8,6))
ax.plot(t, S,'b-', lw=2, label="nicht immune Erwachsene")
ax.plot(t, I,'r--',lw=2, label="Infizierte")
ax.plot(t, R,'g:', lw=2, label="Genesene")
ax.legend(loc="best")
ax.set_xlabel("Zeit/ Tage")
ax.set_ylabel("Erwachsene")
ax.grid(True)
plt.show()       

[OSWALD]

15.3.20.23
Ein wesentlcher Fortschritt in der Erstellung von DGLs ist
die Simulation von DGLs durch Modellbildung, in scipy
minutioes beschrieben.
Hier ein Beispiel , das mein Verständnis
von DGL wesentlich vertieft hat

Eine Modell-Funktion wird in einem 'State Space '
installiert und integriert. Entstanden ist eine
neue DGL mit den Eigenschaften der Modellgleichung.
Ein bedeutender Fortschritt für die Erstellung von
Differentialgleichungen.
Gute Zeit OSWALD

Code: Alles auswählen

import  numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy.integrate import odeint

#tau * dy2/dt2 + 2*zeta*tau*dy/dt + y = Kp*u
#die Modellgleichung

Kp = 10    #  gain
tau = 1    # Zeitkonstante
zeta = 0.25 # Dämpfungsfaktor
theta = 0 # keine Zeitverzögerung,also MULL
du = 1.0    # Verwandlung in u ,u = geeigneterweise glatte Funktion von t 

# (1) Transfer Function
num = [Kp]
den = [tau**2,2*zeta*tau,1]
sys1 = signal.TransferFunction(num,den)
t1,y1 = signal.step(sys1)
 
 

# (2) State Space
A = [[0.0,1.0],[-1.0/tau**2,-2.0*zeta/tau]]  # arraylike
B = [[0.0],[Kp/tau**2]]                      # arraylike 
C = [0.0,0.0]                                # arraylike
D = 0.0                                      # arraylike
 

    
sys2 = signal.StateSpace(A,B,C,D)
#print(sys2)
#print()
t2,y2 = signal.step(sys2)
#print(t2,y2)

#   Ab hier wird die Modellgleichung  von'model3' 
#   transferiert und integriert.   
 
def model3(x,t):
    y = x[0]
    dydt = x[1]
    dy2dt2 = (-2.0*zeta*tau*dydt - y + Kp*du)/tau**2
    return [dydt,dy2dt2] #
 

t3 = np.linspace(0,30,30)
x3 = odeint(model3,[0,0],t3)  #hier wird integriert
y3 = x3[:,0]
 
plt.figure(1)
plt.plot(t1,y1*du,'b-',linewidth=1,label='Transfer Funktion')
 
plt.plot(t2,y2*du,'g:',linewidth=4,label='State Space')
plt.plot(t3,y3,'r*',linewidth=1,label='ODE Integrator')
y_ss = Kp * du
 

plt.plot([0,max(t1)],[y_ss,y_ss],'k:')
plt.xlim([0,max(t1)])
plt.xlabel('Zeit')
plt.ylabel('Antwort(y)')
plt.legend(loc='best')
plt.savefig('2nd_order.png')
plt.grid()
plt.title("Modell Transfer generiert Gfl: Funktionsverläufe identisch ")
plt.show()



[code]

[OSWALD]

18.3.2023
Hier zwei leicht nachvollziehbare Beispiele
für unterschiedliche Lösungen von Dgls.
Gute Zeit Oswald

Code: Alles auswählen

#   Chemie: Reaktionskinetik
import matplotlib.pyplot as plt
import numpy as np
# Lösung mit dem Euler Verfahren
k = 4 # Kinetischer Koeffizient          # Reaktionsgeschwindikeit
Konzentrationen = []; Zeitpunkte = []    # zwei Listen als 'Datenspeicher'
#
#Speicherplatz für Konzentrationen und Zeiten
Konzentration = 1 #Anfangsbedingung       Anfangfskonzentration
tmax = 1 # Zeitdauer der Simulation: 1 Sek   L
dt = 0.01 # Zeitschritt 
Zeit = 0
for i in range(int(tmax/dt)): 
    Konzentrationen.append(Konzentration)
    Zeitpunkte.append(Zeit)
    Konzentration = Konzentration - k * dt * Konzentration # !!!!
    Zeit += dt # !!!!

plt.plot(Zeitpunkte,Konzentrationen)
plt.ylabel('Konzentration [mol/L]')
plt.xlabel('Zeit [s]')
plt.title("Reaktionskinetik nach Euler")
plt.grid()
plt.show()



# Lösung mit einem importieren ODE-Solver aus scipy
tmax = 1 # Zeitdauer der Simulation: in  1 Sek   --> 0
dt = 0.01 # Zeitschritt
k = 4 # Kinetischer Koeffizient         #s.o.
from scipy.integrate import solve_ivp
Stuetzstellen = np.arange(0,tmax,dt)
Konz_ini = 1
def reaktion(t, x): return - k * x # Die Differentialgleichung
### Und hier der zentrale Schritt
sol = solve_ivp(reaktion, [0, tmax], [Konz_ini], t_eval=Stuetzstellen)
### Ende der Rechnung
plt.plot(Stuetzstellen,sol.y[0])
plt.ylabel('Konzentration [mol/L]')
plt.xlabel('Zeit [s]')

plt.title("Reaktionskinetik nach scipy solve_ivp ")
plt.grid()
plt.show()

[OSWALD]

]21.3.2023 Zuletzt noch zwei einfache gekoppelte#DGL-Systeme.
Hierbei geht es um ballistische Funktionen.
gute Zeit Oswald

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 = 3.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, 3))
plt.plot(U[:,0], yi)
plt.plot([u1_H], [H], 'ro')
plt.xlabel('$u_1$')
plt.ylabel('$y$')
plt.show();
 
 
# 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 = 1

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

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

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])
    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[:,2], yi)
plt.plot([u1_H], [H], 'ro')
plt.xlabel('$u_1$')
plt.ylabel('$y$')

 
plt.title("Der rote Punkt ist das Ziel")

plt.show()