from sympy import *
import numpy as np
from matplotlib import pyplot as plt
init_printing()
# Differentialgleichungen höherer Ordnung
x = Symbol('x')
y = Function('y')
dgl = Eq(y(x).diff(x,2), -y(x).diff(x) - y(x))
dgl
lsg = dsolve(dgl)
lsg
mu = S('mu')
char_poly = dgl.subs(y(x), exp(mu*x)).doit()
char_poly
solve(char_poly)
Anfangsbedingungen $y(1)=2$, $y'(1)=3$
lsg = dsolve(dgl, ics={y(1): 2, y(x).diff(x).subs(x, 1): 3})
lsg
phi = lsg.rhs
plot(phi);
nu = S('nu')
dgl = Eq(y(x).diff(x,2) + y(x).diff(x)/x + (1 - nu**2/x**2)*y(x), 0)
dgl
lsg = dsolve(dgl)
lsg
phi = lsg.rhs
print(phi)
C1*besselj(sqrt(nu**2), x) + C2*bessely(sqrt(nu**2), x)
$J_\nu$: Besselfunktion erster Art der Ordnung $\nu$
$Y_\nu$: Besselfunktion zweiter Art der Ordnung $\nu$, Webersche Funktion
nu_n = 2
plot(besselj(nu_n, x), bessely(nu_n, x), (x, 0, 60), ylim=(-3,1));
Wronskische
M = Matrix(2,2,[besselj(nu, x), bessely(nu,x), besselj(nu,x).diff(x), bessely(nu,x).diff(x)])
M
M.det()
W = M.det().simplify()
W
plot(W.subs(nu, 3), (x, .2, 10));
Welche Funktion das ist, sehen wir mit der Reihenentwicklung
Allerdings müssen wir diese Entwicklung selber implementieren, weil es sich nicht um eine Taylorreihe handelt
ordnung = 4
def Bj_ser(nu, x):
k = S('k')
a = (-1)**k * (x/2)**(2*k) / factorial(k) / gamma(nu+1+k)
res = 0
for kk in range(ordnung):
res += a.subs(k, kk)
return res * (x/2)**nu
Bj_ser(nu, x)
Bj_ser(-5.5, 2)
besselj(-5.5, 2).n()
mu = Wild('mu')
z = Wild('z')
connection_formula = {bessely(mu,z): (besselj(mu,z)*cos(mu*pi)-besselj(-mu,x))/sin(mu*pi)}
connection_formula
W
W1 = W.replace(bessely(mu,z), connection_formula[bessely(mu,z)])
W1
W2 = W1.replace(besselj, Bj_ser)
W2
W3 = W2.expand(trig=True) + O(x**ordnung)
W3
simplify(W3.removeO())
Im richtigen Leben müsste man z.B. mit ordnung=8 und ordnung=20 wiederholen
y = Function('y')
w = Function('w')
t = Symbol('t', positive=True)
a = Rational(1,7)
b = Rational(1,6)
dgl1 = Eq(y(t).diff(t,2), a*w(t) - (1+a)*y(t))
dgl2 = Eq(w(t).diff(t,2), b*y(t) - (1+b)*w(t))
dgs = {dgl1, dgl2}
dgs
lsg = dsolve(dgs)
lsg
phi_c = lsg[0].rhs
psi_c = lsg[1].rhs
Bestimme reelles Fundamentalsystem
C = symbols('C1, C2, C3, C4')
ersetzungen = []
for j in range(4):
ers = {C[j]: 1}
for i in range(4):
if j != i:
ers[C[i]] = 0
ersetzungen.append(ers)
ersetzungen
fus = []
for s in ersetzungen:
fus.append([phi_c.subs(s),
psi_c.subs(s)])
fus
f = fus[0][0] + fus[2][0]
g = fus[0][1] + fus[2][1]
f, g
plot(f,g+3, (t,-100,100));