sympy

HM3 Sympy

Aufgabe 1

# Aufgabe 1
import sympy as sp

# Koordinaten in kartesischen Variablen
x, t = sp.symbols('x t', real=True)

# Koeffizienten bestimmen
M = 2*t*x**3 - 2*t
N = 3*t**2*x**2

#M = 3*t*x**2
#N = t*(2*x**3-2)

# Exaktheit prüfen: dM/dx = dN/dt ?
dMdx = sp.diff(M, x)
dNdt = sp.diff(N, t)
is_exact = sp.simplify(dMdx - dNdt) == 0
print("Die DGL ist exakt:", is_exact)

# Potenzialfunktion Phi finden
Phi = sp.integrate(M, t) + sp.integrate(N - sp.diff(sp.integrate(M, t), x), x)
print("Potenzialfunktion Phi:", Phi)
# Lösung der DGL bestimmen
solution = sp.Eq(Phi, sp.Symbol('C'))
print("Allgemeine Lösung:", solution)

# Test
# N
Phi.diff(x)-M
# erwartet: 0
# M
Phi.diff(t)-N
# erwartet: 0

# Anfangswertproblem lösen
C_value = sp.solve(solution.subs([(t, 1), (x, 0)]), sp.Symbol('C'))
final_solution = solution.subs(sp.Symbol('C'), C_value[0])
print("Lösung des Anfangswertproblems:", final_solution)

# Grenzwert für t -> ∞ bestimmen
x_infinity = sp.limit(final_solution.rhs, t, sp.oo)
print("Grenzwert der Lösung für t -> ∞:", x_infinity)

Aufgabe 2

# Aufgabe 2
# System von ODE
import sympy as sp

# variables
x,y,c1,c2 = sp.symbols('x y c_1 c_2')
t = sp.symbols('t', positive = True)
A=sp.Matrix([(3,-1),(5, -1)])
#A=sp.Matrix([(0,-1),(5, 4)])


# calculate only one eigenvalue
# the other one is complex conjugate
EV=(sp.trace(A)+sp.sqrt(sp.trace(A)**2-4*sp.det(A)))/2
print(EV)

# calculate one eigenvector
X=sp.Matrix((x,y))
M=A-EV*sp.eye(2)
M1=M*X

# solve only first row, set x=1
y=sp.solve(M1[0], y)
# choose x=1 and calculate y
X=sp.Matrix((1,(y[0].subs(x,1))))

# express all variables (a=re(X),b=im(X)) and (re(EV), im(EV) in the real solution of the system of ODEs
#a=X.as_real_imag()[0]
#b=X.as_real_imag()[1]
a, b = X.as_real_imag()
print(a)
print(b)
#ev_real=EV.as_real_imag()[0]
#ev_imag=EV.as_real_imag()[1]
ev_real, ev_imag = EV.as_real_imag()
# solution
x1=sp.exp(ev_real*t)*(a*sp.cos(ev_imag*t)-b*sp.sin(ev_imag*t))
x2=sp.exp(ev_real*t)*(a*sp.sin(ev_imag*t)+b*sp.cos(ev_imag*t))
matrix = sp.Matrix([(x1,x2)])
#print(matrix)
solution=c1*x1+c2*x2
#print(solution)
solution

# partikuläre Lösung
B=sp.Matrix((1,0))
xp=-A.inv()*B
xp
print (xp)


# Test ob EV korrekt ist
assert sp.simplify(M*X)==sp.Matrix([0,0]), 'Eigenvektor nok'
matrix

A=sp.Matrix([(3,-1),(5, -1)])
λ1,λ2=A.eigenvals()
print(λ1,λ2)
v1,v2=A.eigenvects()
print(v1,v2)

sol=5/2*v1[2][0]*sp.exp(λ1*t)+5/2*v2[2][0]*sp.exp(λ2*t)
sol1, sol2 = sol.as_real_imag()
sol

# short sympy solution
A=sp.Matrix([(3,-1),(5, -1)])

# Berechne Eigenwerte und Eigenvektoren
eigen_data = A.eigenvects()

# Skaliere die Eigenvektoren so, dass der erste Eintrag 1 wird
# man teilt durch den ersten Eintrag (hier komplex)
eigenvectors = []
for eigenvalue, multiplicity, eigenvecs in eigen_data:
    for vec in eigenvecs:
        if vec[0] != 0:  # Verhindern, dass durch 0 geteilt wird
            vec = vec / vec[0]  # Normierung auf ersten Eintrag = 1
            
        eigenvectors.append((eigenvalue, vec))

# Komplexe Lösung
λ1, v1  = eigenvectors[0]
#λ2, v2  = eigenvectors[1]
v1 = A.eigenvects()[0][2][0]/A.eigenvects()[0][2][0][0]
sol1=v1*sp.exp(λ1*t) # es genügt eine komplexe Lösung zu betrachten; Rotor

# Umschreiben als Matrix von reellen Lösungen
sol_list=sol1.as_real_imag() # ergibt eine Liste von Vektoren
matrix=sp.Matrix([sol_list]) # als Matrix
#sp.simplify(matrix/sp.exp(t))
matrix

# noch kürzer
# Eigenvektor und Eigenwert
v1 = A.eigenvects()[0][2][0]/A.eigenvects()[0][2][0][0] # erster Eintrag auf 1 normiert
λ1, _ = A.eigenvals() 
# reelles System
sp.Matrix([(v1*sp.exp(λ1*t)).as_real_imag()])

Aufgabe 3

# Aufgabe 3
import sympy as sp

# Parabolic PDE

# Define symbols
x, t, k, a = sp.symbols('x t k a', real=True, positive=True)

# Fourier coefficients for f(x) = 1 on [0, pi]
b_k = (2/sp.pi) * sp.integrate(sp.sin(k*x), (x, 0, sp.pi))
b_k = sp.simplify(b_k)

# Solution for u(0,x) = sin(kx)
u_k = sp.exp(-a * k**2 * t) * sp.sin(k*x)

# General solution for u(0,x) = 1 (sum over odd k)
n = sp.symbols('n', integer=True, positive=True)
k_n = 2*n - 1  # Only odd indices
u_series = sp.Sum((4 / (k_n * sp.pi)) * sp.exp(-a * k_n**2 * t) * sp.sin(k_n * x), (n, 1, sp.oo))

# Display results
print("Fourier coefficient b_k:", b_k)
print("Solution for u(0,x) = sin(kx):", u_k)
print("General solution u(t,x):", sp.simplify(u_series))
b_k
u_k
u_series

Aufgabe 4

import sympy as sp

# Definition der Variablen
x, y, z, t, beta = sp.symbols('x y z t beta')

# Definition des Vektorfeldes F
F = sp.Matrix([
    y**2 - 2*beta*z*sp.exp(x*z),
    2*x*y + 3*z,
    3*y - 2*x*sp.exp(x*z)
])

# Berechnung der Rotation von F
rot_F = sp.simplify(sp.Matrix([sp.diff(F[2], y) - sp.diff(F[1], z),
                               sp.diff(F[0], z) - sp.diff(F[2], x),
                               sp.diff(F[1], x) - sp.diff(F[0], y)]))

# Lösung für beta, sodass rot F = 0
beta_solution = sp.solve(rot_F, beta)

# Einsetzen von beta = 1
F_conservative = F.subs(beta, 1)

# Bestimmung des Potentials phi
phi_x = sp.integrate(F_conservative[0], x)
phi = phi_x + sp.Function('g')(y, z)

# Ableitung von phi nach y und Vergleich mit F[1]
g_yz = sp.simplify(sp.diff(phi, y) - F_conservative[1])
g_y_solution = sp.integrate(3*z, y)  # Lösung für g(y, z)
phi = phi_x + g_y_solution + sp.Function('h')(z)

# Ableitung von phi nach z und Vergleich mit F[2]
h_z = sp.simplify(sp.diff(phi, z) - F_conservative[2])
h_solution = sp.integrate(0, z)  # h(z) ist konstant
phi = phi.subs(sp.Function('h')(z), h_solution)

# Definition der Kurve gamma
gamma = sp.Matrix([
    sp.cos(t) - t/(2*sp.pi),
    sp.cos(t),
    sp.sin(t)
])

# Evaluierung des Potentials an den Endpunkten
phi_start = phi.subs({x: gamma[0].subs(t, 0), y: gamma[1].subs(t, 0), z: gamma[2].subs(t, 0)})
phi_end = phi.subs({x: gamma[0].subs(t, 2*sp.pi), y: gamma[1].subs(t, 2*sp.pi), z: gamma[2].subs(t, 2*sp.pi)})

# Berechnung des Kurvenintegrals
integral_value = phi_end - phi_start

# Ausgabe der Ergebnisse
print("Lösungen für beta, sodass F konservativ ist:")
print(beta_solution)
print("\nPotential phi(x,y,z) für beta=1:")
print(phi)
print("\nWert des Kurvenintegrals:")
print(integral_value)

Aufgabe 5

# Aufgabe 5
# mühsam von Hand
# Parametrisierung
u,v = sp.symbols('u v')
#g_u = sp.Matrix([sp.exp(u)-1,0,u])
#g_u = sp.Matrix([u**2,0,u])
#g_u=sp.Matrix([u**2-1,0,u])
g_u=sp.Matrix([u,0,-u**2])
g_uv =sp.Matrix([ g_u[0]*sp.cos(v), g_u[0]*sp.sin (v), u])
f_u=g_uv.diff(u)
f_v=g_uv.diff(v)
n=sp.simplify(f_u.cross(f_v))
print('Normal vector: ')
n

# Substitution F(x,y,z)->F(gamma(u,v))
x,y,z = sp.symbols('x y z')
#F_xyz=sp.Matrix([x+y,y-x,2]) # Typical vector field with div F = 2
#F_xyz=sp.Matrix([x,y,z])
F_xyz=-sp.Matrix([y,-x,z])
F_uv=F_xyz.subs(x,g_uv[0]).subs(y,g_uv[1]).subs(z,g_uv[2])
print('Parametrisiertes Vektorfeld: ')
F_uv

integrand=F_uv.dot(n)
#integrand=F_uv.dot(n).simplify()
print(integrand)
print(integrand.simplify())
integrand.simplify()

sp.simplify(sp.integrate(integrand,(u,0,2), (v,0,2*sp.pi)))

Aufgabe 6

import sympy as sp

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

# Vektorfelder definieren
F = sp.Matrix([2*y*sp.exp(z) + x, y - 2*x*z + 9, z*x + sp.sin(y)])
G = sp.Matrix([-y*z, x*z, z])

# Divergenz von F berechnen
div_F = sp.simplify(sp.diff(F[0], x) + sp.diff(F[1], y) + sp.diff(F[2], z))

# Volumenelement in Zylinderkoordinaten
dV = r * sp.diff(x, x) * sp.diff(r, r) * sp.diff(phi, phi)

# Volumenintegral für den Fluss
flux_integral = sp.integrate(
    sp.integrate(
        sp.integrate((x + 2) * r, (r, 0, x)),
        (phi, 0, 2*sp.pi)
    ),
    (x, 1, 3)
)

# Rotation von G berechnen
rot_G = sp.simplify(sp.Matrix([
    sp.diff(G[2], y) - sp.diff(G[1], z),
    sp.diff(G[0], z) - sp.diff(G[2], x),
    sp.diff(G[1], x) - sp.diff(G[0], y)
]))

# Divergenz der Rotation von G
div_rot_G = sp.simplify(sp.diff(rot_G[0], x) + sp.diff(rot_G[1], y) + sp.diff(rot_G[2], z))

# Ausgabe
print("Divergenz von F:", div_F)
print("Fluss durch K:", flux_integral)
print("Rotation von G:", rot_G)
print("Divergenz der Rotation von G:", div_rot_G)

import sympy as sp

# Parameter für den Kreisrand
r, t = sp.symbols('r t')
r = 2
x = r * sp.cos(t)
y = r * sp.sin(t)
z = 0  # Kreis liegt in der xy-Ebene



# Vektorfeld F
F = sp.Matrix([x - y, x + y, z + x*y])

# Tangentialvektor dr/dt
r_t = sp.Matrix([x, y, z])
dr_dt = sp.diff(r_t, t)

# Skalarprodukt F · dr/dt
integrand = F.dot(dr_dt)

# Kurvenintegral entlang des Kreises
integral = sp.integrate(integrand, (t, sp.pi, 0))

# Kurvenintegral 

# Ergebnis ausgeben
sp.pprint(integral.simplify())