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

[OSWALD]

11.04.2023
Die Interpolation ist ein mathematisches Werkzeug,
das zur Schätzung des Wertes einer Funktion zwischen zwei bekannten Datenpunkten verwendet wird.
Es kann in einer Vielzahl von Kontexten verwendet werden,
z. B. in der Mathematik, Informatik und Statistik.

Die Genauigkeit der Schätzungen hängt von der Qualität der Daten ab,
die für die Schätzungen verwendet werden.
Mit Python und seinen Bibliothen kann die Qualität der Schätzungen
erheblich verbessert werden.

Ich habe im Internet von "Octave" gelesen, das besonders
geeignet sein soll, die Interpolate-Probleme richtig zu lösen.
Das hat mich dazu veranlasst, mich zunächst einmal bei numpy und
scipy zu informieren, was ich dann auch gründlich getan habe.
Wenn ich glaube , diese komplizierte Materie etwa verstanden zu
haben, werde ich Vergleiche mit 'Octave' anstellen.
Jetzt aber erst einmal das, was scipy und numpy in der Interpolation leisten.
Mit ersten Beipielen werde ich In das Thema 'Interpolation' einsteigen.

Gute Zeit OSWALD

[OSWALD]

11.4.23
Hier ein erstes einfachesBeipiel.
Es sollte selbsterklärend sein., dient nur als Illustration
Damit ist aber noch lange nicht gesagt,
ob damit in der Praxis auch ein gutes Ergebnis
erzielt werden kann.
Die Visualisierung dient neben der notwendigen
Mathematik nur als ein möglches Modell.

Code: Alles auswählen

#Univariate    1D Interpolation  
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
x = np.arange(0, 20)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y)  #Interpolation erfolgt über `interp1d`
xnew = np.arange(4.0 , 5., 0.01)

#Achtung :die Interpolation erfolgt hier  von 4.0 bis 5.0 
#        nächster  Abstand  zwischen zwei  Datenpunkte.
ynew = f(xnew)   
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.title("Univariate 1D Interpolation  ") 
plt.grid()
plt.show()

[OSWALD]

In loser Reihenfolge weitere Beispiele für eInzelne
Methoden. Lagrange Methode und Beiepiel für
Suche nach geeigneter Funktion für Interpolation.
OSWALD

Code: Alles auswählen

import numpy as np
from scipy.interpolate import lagrange
x = np.array([0, 1, 2])

y = x**3         # Funktionen stimmen überein
#y = x**3 +3     # Funktion zeigen Parallelverschiebung 
# hier 2 unterschieliche Funktionen  im Versuch

poly = lagrange(x, y)

#nach der Ableitung erhält man   f' =  3*x**2 -2*x                              

 
 
from numpy.polynomial.polynomial import Polynomial
Polynomial(poly.coef[::-1]).coef
#array([ 0., -2.,  3.])
import matplotlib.pyplot as plt
x_new = np.arange(0, 2.1, 0.1)
plt.scatter(x, y, label='Datenpunkte')
plt.plot(x_new, Polynomial(poly.coef[::-1])(x_new), label='Polynomial')
plt.plot(x_new, 3*x_new**2 - 2*x_new + 0*x_new,
         label=r"$3 x^2 - 2 x$", linestyle='-.')
plt.legend()
plt.title("Lagrange- Polynome enthalten immer 3 Stützpunkte ") 
plt.show()

[OSWALD]

12.04.23
Das Taylor-Polynom ist nur eine der von bedeutenden
Mathematikern entwickelten Methoden zur Lösung des Interpolation-Problems.
IN scipy werden sie lückenlos beschrieben.

Gute Zeit OSWALD

Code: Alles auswählen

#Das Taylorpolynom ist eine Näherung für Funktionswerte
# von f = np.sIn(x)  in der Nähe vom Entwicklungspunkt a. 


import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import approximate_taylor_polynomial
x = np.linspace(-10.0, 10.0, num=30)
print(x)   #Stützunkt-Interrvall:   ( hier  -0.344 bis  +0.344)

plt.plot(x, np.sin(x), label="sin-Kurve")
for degree in np.arange(1, 5, step=2):
    sin_taylor = approximate_taylor_polynomial(np.sin, 1, degree, 1,
                                               order=degree + 2)
    plt.plot(x, sin_taylor(x), label=f"degree={degree}")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left',
           borderaxespad=0.0, shadow=True)
plt.tight_layout()
plt.axis([-10, 10, -10, 10])
plt.title(" Taylor-Polynom Näherung von f''' bei Punkt x  ") 
plt.grid()
plt.show()

[OSWALD]

15.04. .2023
ich bin immer noch damit beschäftigt das
Interpolations-Problem wirklich zu verstehen.
Dabei hat mir auch 'octave' nicht geholfen.
'Octave ' scheint mir für eher für Leute geeignet zu sein,
die eher ein fertiges Werkzeug für die Anwendung
in der Praxis suchen. Für Lehrlinge weniger geeignet.

Ich aber bin immer noch damit beschäftigt das
Interpolations-Problem wirklich zu verstehen.

Nach diesem Zitat eines Profis weiß ich endlich, was ein
SPline ist.

'Spline' — Dies bedeutet einfach ein stückweises Polynom
vom Grad k , das k-1 mal kontinuierlich differenzierbar ist

Dafür habe ich hier zwei Beispiele :


code]

import numpy as np
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt
x = np.arange(10)
y = np.sin(x)
cs = CubicSpline(x, y)
xs = np.arange(-0.5, 9.6, 0.01)
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.plot(x, y, 'o', label='Daten')
ax.plot(xs, np.sin(xs), label='true')
ax.plot(xs, cs(xs), label="S")
ax.plot(xs, cs(xs, 1), label="S'")
ax.plot(xs, cs(xs, 2), label="S''")
ax.plot(xs, cs(xs, 3), label="S'''")
ax.set_xlim(0, 10)
ax.legend(loc='lower left', ncol=2)

plt.title("Die verschiedenen Ableitungen von sin(xs) ")
plt.show()

#../../_images/scipy-interpolate-CubicSpline-1_00_00.png
#Im zweiten Beispiel wird der Einheitskreis mit einem Spline interpoliert.
# Es wird eine periodische Randbedingung verwendet.
#Sie können sehen, dass die ersten Ableitungswerte,
# ds/dx=0, ds/dy=1 am periodischen Punkt (1, 0) korrekt berechnet werden.
#Beachten Sie, dass ein Kreis
# nicht exakt durch einen kubischen Spline dargestellt werden kann.
#Um die Genauigkeit zu erhöhen, wären mehr Haltepunkte erforderlich.

theta = 2 * np.pi * np.linspace(0, 1, 5)
y = np.c_[np.cos(theta), np.sin(theta)]
cs = CubicSpline(theta, y, bc_type='periodic')
print("ds/dx={:.1f} ds/dy={:.1f}".format(cs(0, 1)[0], cs(0, 1)[1]))
#cs/dx=0.0
#cs/dy=1.0
s = 2 * np.pi * np.linspace(0, 3, 100)
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.plot(y[:, 0], y[:, 1], 'o', label='Daten')
ax.plot(np.cos(xs), np.sin(xs), label="cos(xs)**2 + sin(xs)**2 =1")
ax.plot(cs(xs)[:, 0], cs(xs)[:, 1], label='Spline')
ax.axes.set_aspect('equal')
ax.legend(loc='center')
plt.title("Der CubicSpline als 'Einheitskreis")
plt.show()


[/code]

[OSWALD]

16.04.2023
Hier nun eine weitere Darstellung
einer cubi Spline interpolation.
und drei Punkten.
Würde man diese Interpolation rein mathematisch
vornehmen , wäre das eine Plagerei.
Deshalb sind die graphischen Darstellungen nicht
nur schöne Bildchen', sondern durchaus sinnvol.
Man kann hier zwar die Koordinaten in etwa ablesen.
Um genaue Koordinaten zu eerhalten , habe ich
jeweils 10 VWerte für x _new und y_ new ausgedruckt.
(Zeile 13).
Gute Zeit OSWALD

Code: Alles auswählen

#ubic spline interpolation
 
from scipy.interpolate import CubicSpline
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('seaborn-poster')
x = [0, 1, 2]
y = [1, 3, 2]

# use bc_type = 'natural' adds the constraints as we described above
f = CubicSpline(x, y, bc_type='natural')
x_new = np.linspace(0, 2 , 10, 0.01)
y_new = f(x_new)
plt.figure(figsize = (8,8))
plt.plot(x_new, y_new, 'b')
plt.plot(x, y, 'ro')
plt.title('Cubic Spline Interpolation')
plt.xlabel('x')
plt.ylabel('y')
   
 
print("X-Koordinaten",x_new)
print("Y-Koordinaten",y_new)
plt.grid()
plt.show()

[OSWALD]

16.04.23
zum Thema intra- extrapolieren ' och eine Bemerkung.
Wenn man in der Zeile 13 den Wert von 10 auf 20 erhöht
wird die Darstellung sichtbar genauergenauer.
Man könnte nun beim Wert 10 eine Extrapolation vornehmen,
inde, man je 2 Koordinaten addiert und dann halbiert.
Das Gleiche erreichen wir, wenn wir die Zahl auf 20 erhöhen.
Mit cubic Spline wird die Interpolation aber vervielfacht, genau,
indem wir unsere Zahl beliebig erhöhen.
Das ist meine aktuelle Vorstellung von 'Interpolation.'.
Ist das korrekt ?

[OSWALD]

17.04.23
zum Themaabschluss noch zwei
interessante Programme,
Ob es für eine erfolgreiche praktische
Anwendung reicht ? Ic glaube es nicht, aber
irgendwann könnte es funktionieren.
Gute Zeit OSWALD

Code: Alles auswählen

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

def find_roots(x, y):
    s = np.abs(np.diff(np.sign(y))).astype(bool)
    return x[:-1][s] + np.diff(x)[s] / (np.abs(y[1:][s] / y[:-1][s]) + 1)

x = [0.1, 1.5, 2.3, 5.5, 6.7, 7]
y = [5, 4, 5, 2, 2, 3]
s = interpolate.UnivariateSpline(x, y, s=0)

xs = np.linspace(min(x), max(x), 500)  # values for x axis
ys = s(xs)

plt.figure()
plt.plot(xs, ys, label='Spline')
for y0 in range(0, 7):
    r = find_roots(xs, ys - y0)
    if len(r) > 0:
        plt.scatter(r, np.repeat(y0, len(r)), label=f'y0 = {y0}')
        for ri in r:
            plt.text(ri, y0, f'  {ri:.3f}', ha='left', va='center')
plt.legend()
plt.xlim(min(x) - 0.2, max(x) + 1)
plt.show()

###############################
#(sollte selbsterklärend sein)
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
x = np.arange(0, 10)
#y = np.exp(-x/3.0)
y  = 2*x**2 - 5*np.cos(x**2) 
f = interpolate.interp1d(x, y, fill_value='extrapolate')
#f = interpolate.interp1d(x, y)

print(f(1))
print(f(3))

[OSWALD]

8.5.2023
Hallo ich bin wieder hier.
Ich habe mich in einem Bereich verlaufen, der
mit 'Mathespielereien nichts mehr zu zun hat.
Gelandet bin ich mitten im Bereich der Statistik, speziell bei
dem Satz von Bayes. Hier habe ich zwar denn Satz von
Bayes als sochen verstanden, aber keineswegs seine
Lösung auf Python-Ebene. Statt dessen hat mich eine
Unmenge von neuen Problemen überrollt.
Es beginnt mit der Auswahl von geeigneten Zufallszahlen,
über Verteilungen und diverse Sätze von Laplace usw usw , etwa
in scipy. Formal habe ich diese alle verstanden , aber wie und wann
ich diese anwenden muss, um richtige Ergebnisse zu erhalten.

An folgendem Beispiel mit Python habe ich experimentiert, um
eventuell zu erkennen, ob die Ergebnisse plausibel oder gar richtig sind.
Begonnen habe ich mit der Erhöhung von 1000 bis zu 1Mio.
Verblüffend fand ich, dass das Programm mit 1 Mio samples
schneller duerchläuft ist , als mit den ursprünglich 1000.
Dann habe ich mir die 'final_data ausdrucken lassen und festgestellt,
dass - nach meinem Verständnis - allenfall die ersten dreI 'Datenpakete'
richtig sein könnten.
Zuletzt kamen die mir gut bekannten 'Cluster. Diese wiederum sind
sichtlich plausibel.
Anbei nun das von mir herunter geladene Beispiel.

Gute Zeit OSWALD

Code: Alles auswählen

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd 
import seaborn as sns
blue_bowl_prob = 0.5 
white_bowl_prob = 0.5 
football_blue = 3 
basket_blue = 0
football_white = 1
basket_white = 4
tot_blue = football_blue+basket_blue
tot_white = football_white+basket_white
p_blue = [football_blue/tot_blue,basket_blue/tot_blue]
p_white = [football_white/tot_white,basket_white/tot_white]
N_array = np.array([10,50,70,100,1000,10**4,10**5,2.5*10**5,5*10**5,7.5*10**5,10**6]).astype(int)
picked_bowl = ['Blue','White']
picked_ball = ['Football','Basketball']
est_prob_array = []
for n in N_array:
    N_events = n
    final_state = []
    for i in range(N_events):
        p_bowl = np.random.choice(picked_bowl,p=[blue_bowl_prob,white_bowl_prob])
        if p_bowl == 'Blue':
            p_ball = np.random.choice(picked_ball,p = p_blue)
        if p_bowl == 'White':
            p_ball = np.random.choice(picked_ball,p = p_white)
        final_state.append(p_bowl+' '+p_ball)
    final_data = pd.DataFrame(final_state).value_counts()
    
    #est_prob = final_data['Blue Football']/(final_data['Blue Football'] 
    #est_prob_array.append(est_prob)
    print(final_data)
 
from sklearn.datasets import make_classification
x,y=make_classification(n_features=2,n_informative=2,n_redundant=0,class_sep=4,n_samples=1000000
                        )
x[:,0],x[:,1] = (x[:,0]*10).astype(int),(x[:,1]*10).astype(int)
plt.figure(figsize=(6,4))
sns.scatterplot(x=x[:,0],y=x[:,1],hue=y,size=1)
plt.xlabel('Feature 1',fontsize=14)
plt.ylabel('Feature 2',fontsize=14)
plt.legend(title='Label',fontsize=10)
plt.show()