# Einzelne Symbole definieren
x = sp.Symbol('x')
y = sp.Symbol('y')
z = sp.Symbol('z')

# Mehrere Symbole auf einmal definieren
a, b, c = sp.symbols('a b c')

# Parameter mit Annahmen definieren (z.B., dass ein Symbol reell oder positiv ist)
t = sp.Symbol('t', real=True)
r = sp.Symbol('r', positive=True)

# Zeit- oder ortsabhängige Funktionen definieren
f = sp.Function('f')(t)
g = sp.Function('g')(x, y)  # Funktion mit zwei Variablen

print("Definierte Symbole:", x, y, z, a, b, c)
print("Symbol mit Annahmen:", r, "ist positiv:", r.is_positive)
print("Funktionssymbole:", f, g)

# Einfache algebraische Ausdrücke
expr1 = x**2 + 2*x + 1
expr2 = y**3 - 3*y + 5

# Trigonometrische Funktionen
expr3 = sin(x)**2 + cos(x)**2
expr4 = sin(2*x) - 2*sin(x)*cos(x)

# Exponential- und Logarithmusfunktionen
expr5 = exp(x) + exp(-x)
expr6 = log(x*y) - log(x) - log(y)

# Ausdrücke anzeigen
print("Quadratischer Ausdruck:", expr1)
print("Trigonometrische Identität:", expr3)
print("Logarithmus-Eigenschaft:", expr6)

# Ausdrücke vereinfachen
simple_expr = simplify(expr3)  # sin²(x) + cos²(x) sollte 1 ergeben
print("Vereinfachung von sin²(x) + cos²(x):", simple_expr)

simple_expr2 = simplify(expr4)  # sin(2x) - 2sin(x)cos(x) sollte 0 ergeben
print("Vereinfachung von sin(2x) - 2sin(x)cos(x):", simple_expr2)

# Ausdrücke ausmultiplizieren
expanded = expand((x + y)**3)
print("Ausmultiplizieren von (x + y)³:", expanded)

# Ausdrücke faktorisieren
factored = factor(x**2 + 2*x + 1)
print("Faktorisieren von x² + 2x + 1:", factored)

# Terme zusammenfassen
collected = sp.collect(x**2 + y*x + x*y + y**2, x)
print("Zusammenfassen nach x in x² + y*x + x*y + y²:", collected)

# Substitution von Werten oder Ausdrücken
expr_with_subs = (x**2 + y).subs(x, 2)
print("(x² + y) mit x=2:", expr_with_subs)

# Komplexere Substitution
expr_subs_expr = (x**2 + y).subs(x, sin(z))
print("(x² + y) mit x=sin(z):", expr_subs_expr)

# Einfacher symbolischer Ausdruck
expr = sin(x)**2 * exp(-x/10)

# Umwandlung in eine Python-Funktion mit lambdify
f_numpy = sp.lambdify(x, expr, "numpy")

# Numerische Auswertung
x_values = np.linspace(0, 10, 100)
y_values = f_numpy(x_values)

# Visualisierung
plt.figure(figsize=(10, 6))
plt.plot(x_values, y_values)
plt.title(rf'$f(x) = {sp.latex(expr)}$')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.show()

# Funktion mit mehreren Variablen
expr_multi = x**2 + 2*y**2 + sin(x*y)
f_multi = sp.lambdify((x, y), expr_multi, "numpy")

# Auswertung für bestimmte Werte
print(f"f(2, 3) = {f_multi(2, 3)}")

# Auswertung über einen Bereich für 3D-Plot
x_vals = np.linspace(-3, 3, 50)
y_vals = np.linspace(-3, 3, 50)
X, Y = np.meshgrid(x_vals, y_vals)
Z = f_multi(X, Y)

# 3D-Visualisierung
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x,y)')
ax.set_title(rf'$f(x,y) = {sp.latex(expr_multi)}$')
plt.show()

# Einfache Gleichung lösen
eq1 = sp.Eq(x**2 - 4, 0)  # x² - 4 = 0
solution1 = sp.solve(eq1, x)
print(f"Lösungen von {eq1}: {solution1}")

# Gleichung direkt als Ausdruck lösen
solution2 = sp.solve(x**3 - 6*x**2 + 11*x - 6, x)
print(f"Lösungen von x³ - 6x² + 11x - 6 = 0: {solution2}")

# Gleichung mit Parametern lösen
k = sp.Symbol('k')
eq_param = sp.Eq(x**2 + k*x + 1, 0)
solution_param = sp.solve(eq_param, x)
print(f"Lösungen von {eq_param}:\n{solution_param}")

# System von Gleichungen lösen
eq_system = [sp.Eq(x + y, 2), sp.Eq(x - y, 4)]
solution_system = sp.solve(eq_system, (x, y))
print(f"Lösung des Gleichungssystems {eq_system}: {solution_system}")

# Transzendente Gleichung lösen
eq_trans = sp.Eq(sp.sin(x), sp.cos(x))
solution_trans = sp.solve(eq_trans, x)
print(f"Lösungen von {eq_trans}: {solution_trans}")

# Gleichung numerisch lösen
from sympy import nsolve
eq_numeric = sp.Eq(sp.sin(x) - x**2, 0)
solution_numeric = nsolve(eq_numeric, x, 0.5)  # Startwert 0.5
print(f"Numerische Lösung von {eq_numeric} in der Nähe von 0.5: {solution_numeric}")

# Gleichung mit Ungleichung
from sympy import symbols, solve_univariate_inequality, sin, pi
x = symbols('x')
solution_ineq = solve_univariate_inequality(sin(x) > 0.5, x)
print(f"Lösung der Ungleichung sin(x) > 0.5: {solution_ineq}")

# Verschiedene Methoden zum Berechnen von Ableitungen
x, y = sp.symbols('x y')

# Methode 1: diff()-Funktion
f1 = x**3 + x**2 - 2*x + 1
df1_dx = sp.diff(f1, x)
print(f"f(x) = {f1}")
print(f"f'(x) = {df1_dx}")

# Methode 2: Diff-Operator (für mathematisch ansprechendere Darstellung)
f2 = sp.sin(x) * sp.exp(x)
df2_dx = sp.Derivative(f2, x)
print(f"Formale Darstellung der Ableitung: {df2_dx}")
print(f"Berechnete Ableitung: {df2_dx.doit()}")

# Höhere Ableitungen
f3 = sp.log(x)
df3_dx2 = sp.diff(f3, x, 2)  # Zweite Ableitung
print(f"f(x) = {f3}")
print(f"f''(x) = {df3_dx2}")

# Komplexere Funktionen
f4 = sp.sin(x**2) / x
df4_dx = sp.diff(f4, x)
print(f"f(x) = {f4}")
print(f"f'(x) = {df4_dx}")
print(f"Vereinfachte f'(x) = {sp.simplify(df4_dx)}")

# Ableitung an einer bestimmten Stelle auswerten
f5 = x**4 - 4*x**3 + 6*x**2
df5_dx = sp.diff(f5, x)
df5_at_2 = df5_dx.subs(x, 2)
print(f"f(x) = {f5}")
print(f"f'(x) = {df5_dx}")
print(f"f'(2) = {df5_at_2}")

# Ableitung implizit definierter Funktionen
# Für die implizite Funktion x²+y²=1
implicit_eq = x**2 + y**2 - 1
# Um dy/dx zu finden, verwenden wir die implizite Differentiation
# Total differential: 2x + 2y·(dy/dx) = 0
# Umstellung: dy/dx = -x/y
dy_dx = -sp.diff(implicit_eq, x) / sp.diff(implicit_eq, y)
print(f"Für die implizite Funktion {implicit_eq} = 0 ist dy/dx = {dy_dx}")

# Ableitungen von parametrischen Funktionen
t = sp.Symbol('t')
x_t = sp.cos(t)
y_t = sp.sin(t)
# Die Steigung ist dy/dx = (dy/dt)/(dx/dt)
dx_dt = sp.diff(x_t, t)
dy_dt = sp.diff(y_t, t)
slope = dy_dt / dx_dt
print(f"Für die parametrische Kurve x(t)={x_t}, y(t)={y_t} ist dy/dx = {slope}")
print(f"Vereinfacht: dy/dx = {sp.simplify(slope)}")

# Funktion zweier Variablen
x, y, z = sp.symbols('x y z')
f = x**2 * y + sp.sin(y) * sp.exp(x)

# Partielle Ableitung nach x
df_dx = sp.diff(f, x)
print(f"f(x,y) = {f}")
print(f"∂f/∂x = {df_dx}")

# Partielle Ableitung nach y
df_dy = sp.diff(f, y)
print(f"∂f/∂y = {df_dy}")

# Mehrfache partielle Ableitungen
d2f_dx2 = sp.diff(f, x, 2)  # Zweite Ableitung nach x
d2f_dxdy = sp.diff(df_dx, y)  # Gemischte zweite Ableitung
print(f"∂²f/∂x² = {d2f_dx2}")
print(f"∂²f/∂x∂y = {d2f_dxdy}")

# Schwarz'scher Satz: Die Reihenfolge der Differentiation ist egal
d2f_dydx = sp.diff(df_dy, x)
print(f"∂²f/∂y∂x = {d2f_dydx}")
print(f"Sind die gemischten Ableitungen gleich? {d2f_dxdy == d2f_dydx}")

# Gradient berechnen
gradient = [sp.diff(f, var) for var in (x, y)]
print(f"∇f = ({gradient[0]}, {gradient[1]})")

# Gradient als Vektor
from sympy import Matrix
grad_vector = Matrix(gradient)
print("Gradient als Vektor:")
grad_vector

# Gradient an einem bestimmten Punkt auswerten
point = {x: 1, y: 2}
grad_at_point = [expr.subs(point) for expr in gradient]
print(f"∇f an Punkt (1,2) = ({grad_at_point[0]}, {grad_at_point[1]})")

# Gradient für eine Funktion von drei Variablen
g = x**2 + y**2 + z**2 + x*y*z
grad_g = [sp.diff(g, var) for var in (x, y, z)]
print(f"g(x,y,z) = {g}")
print(f"∇g = ({grad_g[0]}, {grad_g[1]}, {grad_g[2]})")

# Anwendung: Richtungsableitung berechnen
# Die Richtungsableitung in Richtung eines Vektors v ist ∇f·v/|v|
v = Matrix([1, 1])  # Richtungsvektor
v_norm = sp.sqrt(v.dot(v))  # Normierung
dir_deriv = (Matrix(gradient).dot(v)) / v_norm
print(f"Richtungsableitung von f in Richtung {v}: {dir_deriv}")

# Beispiel für eine Kurvendiskussion
x = sp.Symbol('x')
f = x**4 - 4*x**3 + 4*x**2 - 4

# 1. Definitionsbereich
# Die Funktion ist polynomiell, also ist der Definitionsbereich ℝ
print("Definitionsbereich: ℝ")

# 2. Nullstellen finden
zeros = sp.solve(f, x)
print(f"Nullstellen von f(x) = {f}:")
for i, zero in enumerate(zeros):
    print(f"x_{i+1} = {zero}")

# 3. Erste Ableitung für Extremwerte
f_prime = sp.diff(f, x)
print(f"f'(x) = {f_prime}")

# Kritische Punkte (f'(x) = 0)
critical_points = sp.solve(f_prime, x)
print("Kritische Punkte (f'(x) = 0):")
for i, cp in enumerate(critical_points):
    print(f"x_{i+1} = {cp}")

# 4. Zweite Ableitung für Minimum/Maximum-Test
f_double_prime = sp.diff(f, x, 2)
print(f"f''(x) = {f_double_prime}")

# Test der kritischen Punkte
print("\nUntersuchung der kritischen Punkte:")
for cp in critical_points:
    cp_value = float(cp)
    f_at_cp = f.subs(x, cp)
    f_double_prime_at_cp = f_double_prime.subs(x, cp)
    
    if f_double_prime_at_cp > 0:
        extremum_type = "lokales Minimum"
    elif f_double_prime_at_cp < 0:
        extremum_type = "lokales Maximum"
    else:
        extremum_type = "Wende- oder Sattelpunkt (weiterer Test nötig)"
    
    print(f"x = {cp_value:.6f}: f(x) = {float(f_at_cp):.6f}, f''(x) = {float(f_double_prime_at_cp):.6f} => {extremum_type}")

# 5. Wendepunkte (f''(x) = 0)
inflection_points = sp.solve(f_double_prime, x)
print("\nWendepunkte (f''(x) = 0):")
for ip in inflection_points:
    ip_value = float(ip)
    f_at_ip = f.subs(x, ip)
    f_triple_prime_at_ip = sp.diff(f, x, 3).subs(x, ip)
    
    if f_triple_prime_at_ip != 0:
        print(f"x = {ip_value:.6f}: f(x) = {float(f_at_ip):.6f} => Wendepunkt")
    else:
        print(f"x = {ip_value:.6f}: f(x) = {float(f_at_ip):.6f} => Kein Wendepunkt")

# 6. Verhalten im Unendlichen
print("\nVerhalten für x → ±∞:")
leading_term = sp.Poly(f, x).all_coeffs()[0] * x**sp.Poly(f, x).degree()
print(f"Dominanter Term: {leading_term}")
if leading_term.coeff(x**sp.Poly(f, x).degree()) > 0:
    print("f(x) → +∞ für x → ±∞")
else:
    print("f(x) → +∞ für x → +∞")
    print("f(x) → -∞ für x → -∞")

# 7. Grafische Darstellung
import numpy as np
import matplotlib.pyplot as plt

# Funktion in numpy-Format konvertieren
f_numpy = sp.lambdify(x, f, "numpy")
x_vals = np.linspace(-1, 5, 1000)
y_vals = f_numpy(x_vals)

plt.figure(figsize=(12, 8))
plt.plot(x_vals, y_vals, 'b-', linewidth=2, label=f'f(x) = {sp.latex(f)}')

# Nullstellen markieren
for zero in zeros:
    try:
        zero_val = float(zero)
        plt.plot(zero_val, 0, 'ro', markersize=8)
        plt.annotate(f'Nullstelle\n({zero_val:.2f}, 0)', 
                     xy=(zero_val, 0), 
                     xytext=(zero_val+0.3, -1),
                     arrowprops=dict(arrowstyle="->", color="red"))
    except:
        pass  # Für komplexe Nullstellen

# Extrempunkte markieren
for cp in critical_points:
    try:
        cp_val = float(cp)
        f_at_cp = float(f.subs(x, cp))
        
        if float(f_double_prime.subs(x, cp)) > 0:
            marker_color = 'go'  # Grün für Minima
            point_type = 'Minimum'
        else:
            marker_color = 'mo'  # Magenta für Maxima
            point_type = 'Maximum'
            
        plt.plot(cp_val, f_at_cp, marker_color, markersize=8)
        plt.annotate(f'{point_type}\n({cp_val:.2f}, {f_at_cp:.2f})', 
                     xy=(cp_val, f_at_cp), 
                     xytext=(cp_val-0.5, f_at_cp+2),
                     arrowprops=dict(arrowstyle="->"))
    except:
        pass  # Für komplexe kritische Punkte

# Wendepunkte markieren
for ip in inflection_points:
    try:
        ip_val = float(ip)
        f_at_ip = float(f.subs(x, ip))
        
        if float(sp.diff(f, x, 3).subs(x, ip)) != 0:
            plt.plot(ip_val, f_at_ip, 'ys', markersize=8)  # Gelbes Quadrat für Wendepunkte
            plt.annotate(f'Wendepunkt\n({ip_val:.2f}, {f_at_ip:.2f})', 
                        xy=(ip_val, f_at_ip), 
                        xytext=(ip_val+0.5, f_at_ip-2),
                        arrowprops=dict(arrowstyle="->", color="orange"))
    except:
        pass

plt.grid(True)
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
plt.legend()
plt.title('Kurvendiskussion')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.ylim(-5,1)
plt.show()

# 8. Zusammenfassung der Kurvendiskussion
print("\nZusammenfassung der Kurvendiskussion für f(x) =", f)
print(f"- Definitionsbereich: ℝ")
print(f"- Nullstellen: {zeros}")
print(f"- Extrempunkte: {critical_points}")
print(f"- Wendepunkte: {inflection_points}")

# Unbestimmte Integrale
x = sp.Symbol('x')
from IPython.display import display, Math

# Einfache Integrale
f1 = x**2
F1 = sp.integrate(f1, x)
display(Math(r"\int " + sp.latex(f1) + r" \, dx = " + sp.latex(F1)))
display(Math(r"\text{Überprüfung: } \frac{d}{dx}\left(" + sp.latex(F1) + r"\right) = " + sp.latex(sp.diff(F1, x))))

# Trigonometrische Funktionen
f2 = sp.sin(x)
F2 = sp.integrate(f2, x)
display(Math(r"\int " + sp.latex(f2) + r" \, dx = " + sp.latex(F2)))

# Exponentialfunktionen
f3 = sp.exp(x)
F3 = sp.integrate(f3, x)
display(Math(r"\int " + sp.latex(f3) + r" \, dx = " + sp.latex(F3)))

# Logarithmische Funktionen
f4 = 1/x
F4 = sp.integrate(f4, x)
display(Math(r"\int " + sp.latex(f4) + r" \, dx = " + sp.latex(F4)))

# Komplexere Integrale
f5 = x * sp.exp(x**2)
F5 = sp.integrate(f5, x)
display(Math(r"\int " + sp.latex(f5) + r" \, dx = " + sp.latex(F5)))
display(Math(r"\text{Überprüfung: } \frac{d}{dx}\left(" + sp.latex(F5) + r"\right) = " + sp.latex(sp.simplify(sp.diff(F5, x)))))

f6 = sp.sin(x)**2 * sp.cos(x)
F6 = sp.integrate(f6, x)
display(Math(r"\int " + sp.latex(f6) + r" \, dx = " + sp.latex(F6)))
display(Math(r"\text{Vereinfacht: } " + sp.latex(sp.simplify(F6))))

# Substitutionsmethode (automatisch von SymPy angewendet)
f7 = x / (x**2 + 1)
F7 = sp.integrate(f7, x)
display(Math(r"\int " + sp.latex(f7) + r" \, dx = " + sp.latex(F7)))

# Partielle Integration (automatisch von SymPy angewendet)
f8 = x * sp.sin(x)
F8 = sp.integrate(f8, x)
display(Math(r"\int " + sp.latex(f8) + r" \, dx = " + sp.latex(F8)))
display(Math(r"\text{Überprüfung: } \frac{d}{dx}\left(" + sp.latex(F8) + r"\right) = " + sp.latex(sp.simplify(sp.diff(F8, x)))))

# Bestimmte Integrale
a, b = sp.symbols('a b')

# Einfaches bestimmtes Integral
g1 = x**2
G1 = sp.integrate(g1, (x, 0, 1))
display(Math(r"\int_{0}^{1} " + sp.latex(g1) + r" \, dx = " + sp.latex(G1)))

# Bestimmtes Integral mit symbolischen Grenzen
G1_sym = sp.integrate(g1, (x, a, b))
display(Math(r"\int_{" + sp.latex(a) + r"}^{" + sp.latex(b) + r"} " + sp.latex(g1) + r" \, dx = " + sp.latex(G1_sym)))

# Trigonometrisches bestimmtes Integral
g2 = sp.sin(x)
G2 = sp.integrate(g2, (x, 0, sp.pi))
display(Math(r"\int_{0}^{\pi} " + sp.latex(g2) + r" \, dx = " + sp.latex(G2)))

# Andere bestimmte Integrale
g3 = 1/(1 + x**2)
G3 = sp.integrate(g3, (x, 0, 1))
display(Math(r"\int_{0}^{1} " + sp.latex(g3) + r" \, dx = " + sp.latex(G3)))
print(f"≈ {float(G3)}")

# Uneigentliches Integral (unendliche Grenze)
g4 = sp.exp(-x)
G4 = sp.integrate(g4, (x, 0, sp.oo))  # sp.oo steht für Unendlich
display(Math(r"\int_{0}^{\infty} " + sp.latex(g4) + r" \, dx = " + sp.latex(G4)))

# Gaussche Glockenkurve
g5 = sp.exp(-x**2)
G5 = sp.integrate(g5, (x, -sp.oo, sp.oo))
display(Math(r"\int_{-\infty}^{\infty} " + sp.latex(g5) + r" \, dx = " + sp.latex(G5)))

# Integration mit Parameter
t = sp.Symbol('t')
g6 = sp.exp(-t*x)
G6 = sp.integrate(g6, (x, 0, sp.oo))
display(Math(r"\int_{0}^{\infty} " + sp.latex(g6) + r" \, dx = " + sp.latex(G6) + r" \text{ (gültig für Re(t) > 0)}"))

# Piecewise-Funktionen (stückweise definierte Funktionen)
g7 = sp.Piecewise((1, x <= 0), (sp.exp(-x), x > 0))
G7 = sp.integrate(g7, (x, -1, 1))
display(Math(r"\int_{-1}^{1} " + sp.latex(g7) + r" \, dx = " + sp.latex(G7)))

# Nicht-elementare Integrale
x = sp.Symbol('x')
from IPython.display import display, Math, Markdown

# Fehler-Funktion (Error Function)
f1 = sp.exp(-x**2)
F1 = sp.integrate(f1, x)
display(Math(r"\int " + sp.latex(f1) + r" \, dx = " + sp.latex(F1)))

# Lösungen mit speziellen Funktionen
display(Markdown("**Die Fehlerfunktion erf(x) ist definiert als:**"))
display(Math(r"\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} \, dt"))

# Integrale, die elliptische Integrale verwenden
f2 = 1/sp.sqrt(1 - x**4)
F2 = sp.integrate(f2, x)
display(Math(r"\int " + sp.latex(f2) + r" \, dx = " + sp.latex(F2)))

# Fresnel-Integrale
f3 = sp.sin(x**2)
F3 = sp.integrate(f3, x)
display(Math(r"\int " + sp.latex(f3) + r" \, dx = " + sp.latex(F3)))

# Bessel-Funktionen
f4 = sp.sin(x)/x
F4 = sp.integrate(f4, (x, 0, sp.oo))
display(Math(r"\int_{0}^{\infty} " + sp.latex(f4) + r" \, dx = " + sp.latex(F4)))

# Liouville-Theorem: Manche elementare Funktionen haben keine
# elementaren Stammfunktionen
f5 = sp.exp(x**2)
F5 = sp.integrate(f5, x)
display(Math(r"\int " + sp.latex(f5) + r" \, dx = " + sp.latex(F5)))

# Weitere Beispiele für nicht-elementare Integrale
f6 = sp.sin(x)/x
F6 = sp.integrate(f6, x)
display(Math(r"\int " + sp.latex(f6) + r" \, dx = " + sp.latex(F6)))

f7 = 1/sp.log(x)
F7 = sp.integrate(f7, x)
display(Math(r"\int " + sp.latex(f7) + r" \, dx = " + sp.latex(F7)))

# Numerische Auswertung nicht-elementarer Integrale
from sympy.abc import z
from mpmath import quad, exp, sin, sqrt, pi

# Numerische Integration der Gaußschen Glockenkurve
def gaussian(t):
    return exp(-t**2)

numerical_result = quad(gaussian, [-sp.oo, sp.oo])
exact_result = sqrt(pi)

display(Markdown("**Numerische Integration von:**"))
display(Math(r"\int_{-\infty}^{\infty} e^{-x^2} \, dx"))
print(f"Numerisches Ergebnis: {numerical_result}")
print(f"Exaktes Ergebnis: {exact_result}")
print(f"Unterschied: {abs(numerical_result - exact_result)}")

# Beispiel für die Sinc-Funktion
def sinc(t):
    if t == 0:
        return 1
    else:
        return sin(t)/t

numerical_sinc = quad(sinc, [0, 100])  # Integration bis zu einem großen Wert statt ∞

display(Markdown("**Numerische Integration von:**"))
display(Math(r"\int_{0}^{100} \frac{\sin(x)}{x} \, dx = " + str(numerical_sinc)))

# Analytisch
display(Math(r"\text{Analytisches Ergebnis: } \int_{0}^{\infty} \frac{\sin(x)}{x} \, dx = \frac{\pi}{2} = " + str(pi/2)))

# Visualisierung nicht-elementarer Integrale
import numpy as np
import matplotlib.pyplot as plt

# Fehler-Funktion
x_vals = np.linspace(-3, 3, 1000)
erf_vals = np.array([float(sp.erf(val)) for val in x_vals])

plt.figure(figsize=(10, 6))
plt.plot(x_vals, erf_vals, 'b-', linewidth=2)
plt.title('Fehlerfunktion erf(x)')
plt.xlabel('x')
plt.ylabel('erf(x)')
plt.grid(True)
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
plt.show()

# Fresnel-Integral
x_vals = np.linspace(0, 10, 1000)
sin_x2_vals = np.sin(x_vals**2)
# Approximation des Integrals durch kumulative Summe
dx = x_vals[1] - x_vals[0]
fresnel_approx = np.cumsum(sin_x2_vals) * dx

plt.figure(figsize=(10, 6))
plt.plot(x_vals, sin_x2_vals, 'b-', linewidth=1, alpha=0.7, label='sin(x²)')
plt.plot(x_vals, fresnel_approx, 'r-', linewidth=2, label='∫sin(x²)dx (approximiert)')
plt.title('Fresnel-Integral S(x) = ∫₀ˣ sin(t²)dt')
plt.xlabel('x')
plt.grid(True)
plt.legend()
plt.show()

# Mehrfachintegrale mit SymPy - Drei Grundlegende Beispiele
import sympy as sp
from IPython.display import display, Math

# Symbole definieren
x, y, z, r, theta, phi = sp.symbols('x y z r theta phi')

# 1. Einfaches Doppelintegral über ein Rechteck
# ∫₀¹∫₀² x·y dy dx
f1 = x * y
I1 = sp.integrate(f1, (y, 0, 2), (x, 0, 1))
display(Math(r"\int_0^1 \int_0^2 " + sp.latex(f1) + r" \, dy \, dx = " + sp.latex(I1)))

# 2. Doppelintegral in Polarkoordinaten (Fläche einer Kreisscheibe)
# Transformation: x = r·cos(θ), y = r·sin(θ)
# Jacobi-Determinante: r
# ∫₀²ᵖ∫₀¹ r dr dθ
f2 = 1
I2 = sp.integrate(f2 * r, (r, 0, 1), (theta, 0, 2*sp.pi))
display(Math(r"\int_0^{2\pi} \int_0^1 " + sp.latex(f2) + r" \cdot " + sp.latex(r) + r" \, dr \, d\theta = " + sp.latex(I2)))

# 3. Dreifachintegral in Kugelkoordinaten (Volumen einer Kugel)
# Transformation: x = r·sin(φ)·cos(θ), y = r·sin(φ)·sin(θ), z = r·cos(φ)
# Jacobi-Determinante: r²·sin(φ)
# ∫₀²ᵖ∫₀ᵖ∫₀ᴿ r²·sin(φ) dr dφ dθ

R = sp.Symbol('R', positive=True)  # Radius der Kugel
f3 = 1
I3 = sp.integrate(f3 * r**2 * sp.sin(phi), (r, 0, R), (phi, 0, sp.pi), (theta, 0, 2*sp.pi))
display(Math(r"\int_0^{2\pi} \int_0^{\pi} \int_0^R " + sp.latex(f3) + r" \cdot " + sp.latex(r**2 * sp.sin(phi)) + 
         r" \, dr \, d\phi \, d\theta = " + sp.latex(I3)))

# Vergleich mit der bekannten Formel V = (4/3)πR³
expected = (4/3) * sp.pi * R**3
display(Math(r"\text{Vergleich mit } \frac{4}{3}\pi R^3 = " + sp.latex(expected) + r": " + 
         sp.latex(sp.simplify(I3 - expected))))

import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative

# Testfunktion definieren
def f(x):
    return np.sin(x) * np.exp(-0.1 * x)

# Exakte Ableitung (zum Vergleich)
def df_exact(x):
    return np.cos(x) * np.exp(-0.1 * x) - 0.1 * np.sin(x) * np.exp(-0.1 * x)

# 1. Vorwärtsdifferenz
def forward_diff(f, x, h=1e-5):
    return (f(x + h) - f(x)) / h

# 2. Rückwärtsdifferenz
def backward_diff(f, x, h=1e-5):
    return (f(x) - f(x - h)) / h

# 3. Zentrale Differenz
def central_diff(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

# 4. Fünf-Punkte-Formel (höhere Genauigkeit)
def five_point_diff(f, x, h=1e-3):
    return (-f(x + 2*h) + 8*f(x + h) - 8*f(x - h) + f(x - 2*h)) / (12 * h)

# Berechnungen für einen einzelnen Punkt
x0 = 2.0
h_values = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9, 1e-10]

# Exakter Wert
exact_derivative = df_exact(x0)
print(f"Exakte Ableitung bei x = {x0}: {exact_derivative}")

# Fehler für verschiedene Methoden und Schrittweiten berechnen
errors_forward = []
errors_backward = []
errors_central = []
errors_five_point = []

for h in h_values:
    # Numerische Ableitungen berechnen
    fd = forward_diff(f, x0, h)
    bd = backward_diff(f, x0, h)
    cd = central_diff(f, x0, h)
    fpd = five_point_diff(f, x0, h)
    
    # Absoluten Fehler berechnen
    errors_forward.append(abs(fd - exact_derivative))
    errors_backward.append(abs(bd - exact_derivative))
    errors_central.append(abs(cd - exact_derivative))
    errors_five_point.append(abs(fpd - exact_derivative))

# Vergleich der Fehler für verschiedene Schrittweiten visualisieren
plt.figure(figsize=(12, 8))
plt.loglog(h_values, errors_forward, 'o-', label='Vorwärtsdifferenz')
plt.loglog(h_values, errors_backward, 's-', label='Rückwärtsdifferenz')
plt.loglog(h_values, errors_central, '^-', label='Zentrale Differenz')
plt.loglog(h_values, errors_five_point, 'D-', label='Fünf-Punkte-Formel')
plt.grid(True, which="both", ls="--")
plt.xlabel('Schrittweite h')
plt.ylabel('Absoluter Fehler')
plt.title('Fehler der numerischen Differentiationsmethoden in Abhängigkeit von der Schrittweite')
plt.legend()
plt.show()

# Visualisierung der Funktion und ihrer Ableitung
x_values = np.linspace(0, 10, 1000)
f_values = f(x_values)
df_values = df_exact(x_values)

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(x_values, f_values)
plt.grid(True)
plt.title('Funktion f(x) = sin(x) · e^(-0.1x)')
plt.xlabel('x')
plt.ylabel('f(x)')

plt.subplot(1, 2, 2)
plt.plot(x_values, df_values)
plt.grid(True)
plt.title('Ableitung f\'(x)')
plt.xlabel('x')
plt.ylabel('f\'(x)')
plt.tight_layout()
plt.show()

# Demonstration der numerischen Ableitung entlang der gesamten Funktion
h_demo = 1e-4  # Eine für Demonstrationszwecke geeignete Schrittweite
central_diff_values = [central_diff(f, x, h_demo) for x in x_values]
five_point_diff_values = [five_point_diff(f, x, h_demo) for x in x_values]

plt.figure(figsize=(12, 6))
plt.plot(x_values, df_values, 'k-', label='Exakte Ableitung')
plt.plot(x_values, central_diff_values, 'r--', label=f'Zentrale Differenz (h={h_demo})')
plt.plot(x_values, five_point_diff_values, 'g:', label=f'Fünf-Punkte-Formel (h={h_demo})')
plt.grid(True)
plt.title('Vergleich von exakter und numerischer Ableitung')
plt.xlabel('x')
plt.ylabel('f\'(x)')
plt.legend()
plt.show()

# Höhere Ableitungen
def second_derivative(f, x, h=1e-4):
    """Berechnet die zweite Ableitung mit der zentralen Differenzenmethode."""
    return (f(x + h) - 2 * f(x) + f(x - h)) / (h**2)

# Exakte zweite Ableitung
def d2f_exact(x):
    return (-0.1*np.cos(x) - 0.1*np.cos(x) - np.sin(x) - 0.1*(-0.1*np.sin(x))) * np.exp(-0.1*x)

# Vergleich der exakten und numerischen zweiten Ableitung
d2f_values = d2f_exact(x_values)
d2f_numerical = [second_derivative(f, x) for x in x_values]

plt.figure(figsize=(10, 6))
plt.plot(x_values, d2f_values, 'k-', label='Exakte zweite Ableitung')
plt.plot(x_values, d2f_numerical, 'r--', label='Numerische zweite Ableitung')
plt.grid(True)
plt.title('Vergleich von exakter und numerischer zweiter Ableitung')
plt.xlabel('x')
plt.ylabel('f\'\'(x)')
plt.legend()
plt.show()

# Demonstration der numerischen Instabilität bei der Differentiation von verrauschten Daten
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
np.random.seed(42)  # Für Reproduzierbarkeit

# Basisfunktion: f(x) = x²
def f_exact(x):
    return x**2

# Exakte Ableitung: f'(x) = 2x
def df_exact(x):
    return 2*x

# Verrauschte Funktion erstellen
def f_noisy(x, noise_level=0.05):
    if np.isscalar(x):
        return f_exact(x) + noise_level * np.random.randn()
    else:
        return f_exact(x) + noise_level * np.random.randn(len(x))

# Numerische Ableitung mit zentraler Differenz
def numerical_derivative(x_grid, h, noise_level=0.05):
    """
    Berechnet die numerische Ableitung für ein gegebenes Gitter, indem zusätzliche 
    Punkte im Abstand h für die Differenzenbildung verwendet werden.
    
    Args:
        x_grid: Das ursprüngliche Gitter
        h: Die Schrittweite für die Differenzenbildung
        noise_level: Stärke des Rauschens
        
    Returns:
        f_grid: Funktionswerte auf dem Originalgitter (verrauscht)
        df_numerical: Numerische Ableitung auf dem Originalgitter
    """
    n = len(x_grid)
    df_numerical = np.zeros_like(x_grid)
    
    # Erzeuge verrauschte Funktionswerte für das Originalgitter
    np.random.seed(42)  # Gleicher Rauschpegel für alle Durchläufe
    f_grid = f_noisy(x_grid, noise_level)
    
    # Für jeden Punkt im Gitter
    for i in range(n):
        x = x_grid[i]
        
        # Erzeuge zusätzliche Punkte links und rechts im Abstand h
        x_minus = x - h
        x_plus = x + h
        
        # Erzeuge verrauschte Funktionswerte für diese zusätzlichen Punkte
        # Wir verwenden hierfür neue Zufallswerte für das Rauschen
        f_minus = f_noisy(x_minus, noise_level)
        f_plus = f_noisy(x_plus, noise_level)
        
        # Berechne die zentrale Differenz
        df_numerical[i] = (f_plus - f_minus) / (2 * h)
    
    return f_grid, df_numerical

# Bereich und Datenpunkte definieren
x = np.linspace(-5, 5, 200)
f_clean = f_exact(x)
df_clean = df_exact(x)

# Schrittweiten für den Vergleich
step_sizes = [1.0, 0.1, 0.01, 0.001]

# Erstellen der Figur mit 4 Subplots
fig = plt.figure(figsize=(14, 10))
gs = GridSpec(2, 2, figure=fig)

# Ergebnisse für MSE-Berechnung speichern
error_info = []

for i, h in enumerate(step_sizes):
    row, col = divmod(i, 2)
    ax = fig.add_subplot(gs[row, col])
    
    # Berechnung der numerischen Ableitung mit aktueller Schrittweite
    f_with_noise, df_numerical = numerical_derivative(x, h)
    
    # Berechnung des MSE
    mse = np.mean((df_numerical - df_clean)**2)
    error_info.append(f'h = {h}: MSE = {mse:.6f}')
    
    # Original verrauschte Funktion plotten
    ax.plot(x, f_with_noise, 'gray', alpha=0.5, linewidth=1, 
            label='Verrauschte Funktion $f(x) = x^2 + \epsilon$')
    
    # Exakte Ableitung plotten
    ax.plot(x, df_clean, 'b-', label='Exakte Ableitung $f\'(x) = 2x$')
    
    # Numerische Ableitung plotten
    ax.plot(x, df_numerical, 'r-', linewidth=1, alpha=0.8, 
            label=f'Numerische Ableitung (h = {h})')
    
    # Formatierung
    ax.set_title(f'Schrittweite h = {h}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid(True, linestyle='--', alpha=0.6)
    
    # Legende einfügen, wenn genug Platz
    if i == 0:
        ax.legend(loc='upper left')
    
    # Achsenbereiche anpassen
    ax.set_ylim(-12, 12)

# Gesamtbeschriftung
plt.suptitle('Numerische Ableitung von verrauschten Daten: Einfluss der Schrittweite h', fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Platz für den Titel oben

# Ausgabe der Fehlerinformationen
print("Mittlerer quadratischer Fehler (MSE) für verschiedene Schrittweiten:")
for info in error_info:
    print(info)

plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 1. Verwendung von numpy.gradient für Datensätze
# Generieren von Beispieldaten
x = np.linspace(0, 10, 100)
y = np.sin(x) * np.exp(-0.2 * x)

# Berechnung der Ableitung mit numpy.gradient
dy_dx = np.gradient(y, x)

plt.figure(figsize=(10, 6))
plt.plot(x, y, 'b-', label='f(x) = sin(x) · e^(-0.2x)')
plt.plot(x, dy_dx, 'r-', label='df/dx (numerisch)')
plt.grid(True)
plt.legend()
plt.title('Numerische Differentiation mit numpy.gradient')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# 2. Implementierung eigener numerischer Ableitungsfunktionen
def f(x):
    return np.sin(x) * np.exp(-0.2 * x)

# Exakte Ableitung zur Überprüfung
def df_exact(x):
    return np.cos(x) * np.exp(-0.2 * x) - 0.2 * np.sin(x) * np.exp(-0.2 * x)

# Eigene Implementation der zentralen Differenz
def central_diff(f, x, h=1e-6):
    return (f(x + h) - f(x - h)) / (2 * h)

# Berechnung der Ableitung an verschiedenen Punkten
x_points = np.linspace(0, 10, 20)
df_numpy = [central_diff(f, x_i) for x_i in x_points]
df_exact_values = [df_exact(x_i) for x_i in x_points]

plt.figure(figsize=(10, 6))
plt.plot(x, [df_exact(x_i) for x_i in x], 'b-', label='Exakte Ableitung')
plt.plot(x_points, df_numpy, 'ro', label='Zentrale Differenz (NumPy)')
plt.grid(True)
plt.legend()
plt.title('Vergleich von exakter und numerischer Ableitung mit zentraler Differenz')
plt.xlabel('x')
plt.ylabel('df/dx')
plt.show()

# Fehler berechnen
errors = np.abs(np.array(df_numpy) - np.array(df_exact_values))
print(f"Maximaler Fehler mit zentraler Differenz: {np.max(errors):.2e}")
print(f"Mittlerer Fehler mit zentraler Differenz: {np.mean(errors):.2e}")

# 3. Höhere Ableitungen mit NumPy
# Eigene Implementation der zentralen Differenz für zweite Ableitung
def second_derivative(f, x, h=1e-4):
    return (f(x + h) - 2 * f(x) + f(x - h)) / (h**2)

# Berechne zweite Ableitung
d2f_numpy = [second_derivative(f, x_i) for x_i in x_points]

# Exakte zweite Ableitung
def d2f_exact(x):
    return -0.2 * np.cos(x) * np.exp(-0.2 * x) - np.sin(x) * np.exp(-0.2 * x) - 0.2 * (np.cos(x) * np.exp(-0.2 * x) - 0.2 * np.sin(x) * np.exp(-0.2 * x))

d2f_exact_values = [d2f_exact(x_i) for x_i in x_points]
d2f_exact_curve = [d2f_exact(x_i) for x_i in x]

plt.figure(figsize=(10, 6))
plt.plot(x, d2f_exact_curve, 'b-', label='Exakte 2. Ableitung')
plt.plot(x_points, d2f_numpy, 'ro', label='Zweite Ableitung (NumPy)')
plt.grid(True)
plt.legend()
plt.title('Zweite Ableitung mit zentraler Differenz')
plt.xlabel('x')
plt.ylabel('d²f/dx²')
plt.show()

# 4. Numerische Differentiation für mehrdimensionale Daten
# Erstellen eines 2D-Beispieldatensatzes
x, y = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
z = np.sin(np.sqrt(x**2 + y**2))

# Berechnung der Gradienten mit numpy.gradient
dx, dy = np.gradient(z, x[0, :], y[:, 0])

# Visualisierung der Funktion und des Gradienten
plt.figure(figsize=(16, 6))

plt.subplot(1, 3, 1)
plt.contourf(x, y, z, 50, cmap='viridis')
plt.colorbar(label='z = sin(√(x² + y²))')
plt.title('2D-Funktion')
plt.xlabel('x')
plt.ylabel('y')

plt.subplot(1, 3, 2)
plt.contourf(x, y, dx, 50, cmap='viridis')
plt.colorbar(label='∂z/∂x')
plt.title('Partielle Ableitung nach x')
plt.xlabel('x')
plt.ylabel('y')

plt.subplot(1, 3, 3)
plt.contourf(x, y, dy, 50, cmap='viridis')
plt.colorbar(label='∂z/∂y')
plt.title('Partielle Ableitung nach y')
plt.xlabel('x')
plt.ylabel('y')

plt.tight_layout()
plt.show()

# Gradientenfeld visualisieren
plt.figure(figsize=(10, 8))
# Für übersichtlichkeit nur jeden 5. Punkt darstellen
skip = 5
plt.contourf(x, y, z, 20, cmap='viridis', alpha=0.3)
plt.quiver(x[::skip, ::skip], y[::skip, ::skip], 
           dx[::skip, ::skip], dy[::skip, ::skip], 
           color='black', scale=50)
plt.title('Gradientenfeld der Funktion z = sin(√(x² + y²))')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from IPython.display import display, Math

# Funktion definieren, die wir integrieren wollen
def f(x):
    return np.sin(x)**2 * np.exp(-0.1 * x)

# Trapezmethode implementieren
def trapezoid_rule(f, a, b, n):
    """
    Numerische Integration mit der Trapezmethode
    
    Parameter:
    f: Zu integrierende Funktion
    a, b: Integrationsgrenzen
    n: Anzahl der Teilintervalle
    
    Rückgabe:
    Näherungswert des Integrals
    """
    h = (b - a) / n  # Schrittweite
    x = np.linspace(a, b, n+1)  # n+1 Punkte für n Intervalle
    y = f(x)
    
    # Trapezformel: h * [0.5*y₀ + y₁ + y₂ + ... + y_{n-1} + 0.5*y_{n}]
    return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])

# Exakten Wert mit SymPy berechnen
x_sym = sp.Symbol('x')
f_sym = sp.sin(x_sym)**2 * sp.exp(-0.1 * x_sym)
a, b = 0, 5  # Integrationsgrenzen
exact_integral = float(sp.integrate(f_sym, (x_sym, 0, 5)))

# Trapezregel mit verschiedenen Anzahlen von Teilintervallen anwenden
n_intervals = [4, 10, 20, 50]
results = []

for n in n_intervals:
    result = trapezoid_rule(f, a, b, n)
    error = abs(result - exact_integral)
    results.append((n, result, error))
    print(f"n = {n:2d}: Trapezregel = {result:.8f}, Fehler = {error:.8f}")

print(f"Exakter Wert:  {exact_integral:.8f}")

# Visualisierung
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Die Funktion und Trapeze für n=4 darstellen
n_vis = 4
x = np.linspace(a, b, n_vis+1)
y = f(x)
x_fine = np.linspace(a, b, 500)
y_fine = f(x_fine)

# Funktion plotten
ax1.plot(x_fine, y_fine, 'b-', linewidth=2, label='f(x)')
ax1.fill_between(x_fine, 0, y_fine, alpha=0.1, color='blue')
ax1.set_title(f'Trapezmethode mit n = {n_vis} Intervallen')
ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.grid(True)

# Trapeze zeichnen
for i in range(n_vis):
    x_trap = [x[i], x[i], x[i+1], x[i+1], x[i]]
    y_trap = [0, y[i], y[i+1], 0, 0]
    ax1.fill(x_trap, y_trap, 'r', alpha=0.2)
    
# Punkte der Auswertung markieren
ax1.plot(x, y, 'ro', markersize=5)

# Konvergenz der Methode zeigen
ax2.plot([r[0] for r in results], [r[2] for r in results], 'b-o', linewidth=2)
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.set_title('Konvergenz der Trapezmethode')
ax2.set_xlabel('Anzahl der Intervalle n')
ax2.set_ylabel('Fehler')
ax2.grid(True, which="both", ls="--")

# Theoretische Konvergenzrate O(h²) = O(1/n²) anzeigen
n_ref = np.array([min(n_intervals), max(n_intervals)])
err_init = results[0][2] * (n_intervals[0]**2)
err_ref = err_init / (n_ref**2)
ax2.plot(n_ref, err_ref, 'k--', label='O(1/n²)')
ax2.legend()

plt.tight_layout()
plt.show()

# Sympy-Darstellung der Trapezformel
from IPython.display import display, Math
display(Math(r"\int_a^b f(x) \, dx \approx \frac{b-a}{n} \left[ \frac{f(a)}{2} + \sum_{i=1}^{n-1} f(x_i) + \frac{f(b)}{2} \right]"))