In [1]:
from sympy import *
import numpy as np
from matplotlib import pyplot as plt
init_printing()
In [2]:
# Differentialgleichungen höherer Ordnung
In [3]:
x = Symbol('x')
y = Function('y')
In [4]:
dgl = Eq(y(x).diff(x,2), -y(x).diff(x) - y(x))
dgl
Out[4]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = - y{\left(x \right)} - \frac{d}{d x} y{\left(x \right)}$
In [5]:
lsg = dsolve(dgl)
lsg
Out[5]:
$\displaystyle y{\left(x \right)} = \left(C_{1} \sin{\left(\frac{\sqrt{3} x}{2} \right)} + C_{2} \cos{\left(\frac{\sqrt{3} x}{2} \right)}\right) e^{- \frac{x}{2}}$
In [6]:
mu = S('mu')
char_poly = dgl.subs(y(x), exp(mu*x)).doit()
char_poly
Out[6]:
$\displaystyle \mu^{2} e^{\mu x} = - \mu e^{\mu x} - e^{\mu x}$
In [7]:
solve(char_poly)
Out[7]:
$\displaystyle \left[ \left\{ \mu : - \frac{1}{2} - \frac{\sqrt{3} i}{2}\right\}, \ \left\{ \mu : - \frac{1}{2} + \frac{\sqrt{3} i}{2}\right\}\right]$

Anfangsbedingungen $y(1)=2$, $y'(1)=3$

In [8]:
lsg = dsolve(dgl, ics={y(1): 2, y(x).diff(x).subs(x, 1): 3})
lsg
Out[8]:
$\displaystyle y{\left(x \right)} = \left(\left(\frac{6 e^{\frac{1}{2}} \sin{\left(\frac{\sqrt{3}}{2} \right)}}{3 \cos^{2}{\left(\frac{\sqrt{3}}{2} \right)} + 3 \sin^{2}{\left(\frac{\sqrt{3}}{2} \right)}} + \frac{8 \sqrt{3} e^{\frac{1}{2}} \cos{\left(\frac{\sqrt{3}}{2} \right)}}{3 \cos^{2}{\left(\frac{\sqrt{3}}{2} \right)} + 3 \sin^{2}{\left(\frac{\sqrt{3}}{2} \right)}}\right) \sin{\left(\frac{\sqrt{3} x}{2} \right)} + \left(- \frac{8 \sqrt{3} e^{\frac{1}{2}} \sin{\left(\frac{\sqrt{3}}{2} \right)}}{3 \cos^{2}{\left(\frac{\sqrt{3}}{2} \right)} + 3 \sin^{2}{\left(\frac{\sqrt{3}}{2} \right)}} + \frac{6 e^{\frac{1}{2}} \cos{\left(\frac{\sqrt{3}}{2} \right)}}{3 \cos^{2}{\left(\frac{\sqrt{3}}{2} \right)} + 3 \sin^{2}{\left(\frac{\sqrt{3}}{2} \right)}}\right) \cos{\left(\frac{\sqrt{3} x}{2} \right)}\right) e^{- \frac{x}{2}}$
In [9]:
phi = lsg.rhs
plot(phi);

Besselfunktionen¶

In [10]:
nu = S('nu')
In [11]:
dgl = Eq(y(x).diff(x,2) + y(x).diff(x)/x + (1 - nu**2/x**2)*y(x), 0)
dgl
Out[11]:
$\displaystyle \left(- \frac{\nu^{2}}{x^{2}} + 1\right) y{\left(x \right)} + \frac{d^{2}}{d x^{2}} y{\left(x \right)} + \frac{\frac{d}{d x} y{\left(x \right)}}{x} = 0$
In [12]:
lsg = dsolve(dgl)
lsg
Out[12]:
$\displaystyle y{\left(x \right)} = C_{1} J_{\sqrt{\nu^{2}}}\left(x\right) + C_{2} Y_{\sqrt{\nu^{2}}}\left(x\right)$
In [13]:
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

In [14]:
nu_n = 2
plot(besselj(nu_n, x), bessely(nu_n, x), (x, 0, 60), ylim=(-3,1));

Wronskische

In [15]:
M = Matrix(2,2,[besselj(nu, x), bessely(nu,x), besselj(nu,x).diff(x), bessely(nu,x).diff(x)])
M
Out[15]:
$\displaystyle \left[\begin{matrix}J_{\nu}\left(x\right) & Y_{\nu}\left(x\right)\\\frac{J_{\nu - 1}\left(x\right)}{2} - \frac{J_{\nu + 1}\left(x\right)}{2} & \frac{Y_{\nu - 1}\left(x\right)}{2} - \frac{Y_{\nu + 1}\left(x\right)}{2}\end{matrix}\right]$
In [16]:
M.det()
Out[16]:
$\displaystyle \frac{J_{\nu}\left(x\right) Y_{\nu - 1}\left(x\right)}{2} - \frac{J_{\nu}\left(x\right) Y_{\nu + 1}\left(x\right)}{2} - \frac{J_{\nu - 1}\left(x\right) Y_{\nu}\left(x\right)}{2} + \frac{J_{\nu + 1}\left(x\right) Y_{\nu}\left(x\right)}{2}$
In [17]:
W = M.det().simplify()
W
Out[17]:
$\displaystyle - J_{\nu}\left(x\right) Y_{\nu + 1}\left(x\right) + J_{\nu + 1}\left(x\right) Y_{\nu}\left(x\right)$
In [18]:
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

https://dlmf.nist.gov/10.2#E2

Pattern matching¶

In [19]:
ordnung = 4
In [20]:
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
In [21]:
Bj_ser(nu, x)
Out[21]:
$\displaystyle \left(\frac{x}{2}\right)^{\nu} \left(- \frac{x^{6}}{384 \Gamma\left(\nu + 4\right)} + \frac{x^{4}}{32 \Gamma\left(\nu + 3\right)} - \frac{x^{2}}{4 \Gamma\left(\nu + 2\right)} + \frac{1}{\Gamma\left(\nu + 1\right)}\right)$
In [22]:
Bj_ser(-5.5, 2)
Out[22]:
$\displaystyle -20.9631692136963$
In [23]:
besselj(-5.5, 2).n()
Out[23]:
$\displaystyle -20.9781995753435$
In [24]:
mu = Wild('mu')
z = Wild('z')
In [25]:
connection_formula =  {bessely(mu,z): (besselj(mu,z)*cos(mu*pi)-besselj(-mu,x))/sin(mu*pi)}
connection_formula
Out[25]:
$\displaystyle \left\{ Y_{\mu}\left(z\right) : \frac{\cos{\left(\pi \mu \right)} J_{\mu}\left(z\right) - J_{- \mu}\left(x\right)}{\sin{\left(\pi \mu \right)}}\right\}$
In [26]:
W
Out[26]:
$\displaystyle - J_{\nu}\left(x\right) Y_{\nu + 1}\left(x\right) + J_{\nu + 1}\left(x\right) Y_{\nu}\left(x\right)$
In [27]:
W1 = W.replace(bessely(mu,z), connection_formula[bessely(mu,z)])
W1
Out[27]:
$\displaystyle \frac{\left(\cos{\left(\pi \nu \right)} J_{\nu}\left(x\right) - J_{- \nu}\left(x\right)\right) J_{\nu + 1}\left(x\right)}{\sin{\left(\pi \nu \right)}} - \frac{\left(\cos{\left(\pi \left(\nu + 1\right) \right)} J_{\nu + 1}\left(x\right) - J_{- \nu - 1}\left(x\right)\right) J_{\nu}\left(x\right)}{\sin{\left(\pi \left(\nu + 1\right) \right)}}$
In [28]:
W2 = W1.replace(besselj, Bj_ser)
W2
Out[28]:
$\displaystyle - \frac{\left(\frac{x}{2}\right)^{\nu} \left(- \left(\frac{x}{2}\right)^{- \nu - 1} \left(- \frac{x^{6}}{384 \Gamma\left(3 - \nu\right)} + \frac{x^{4}}{32 \Gamma\left(2 - \nu\right)} - \frac{x^{2}}{4 \Gamma\left(1 - \nu\right)} + \frac{1}{\Gamma\left(- \nu\right)}\right) + \left(\frac{x}{2}\right)^{\nu + 1} \left(- \frac{x^{6}}{384 \Gamma\left(\nu + 5\right)} + \frac{x^{4}}{32 \Gamma\left(\nu + 4\right)} - \frac{x^{2}}{4 \Gamma\left(\nu + 3\right)} + \frac{1}{\Gamma\left(\nu + 2\right)}\right) \cos{\left(\pi \left(\nu + 1\right) \right)}\right) \left(- \frac{x^{6}}{384 \Gamma\left(\nu + 4\right)} + \frac{x^{4}}{32 \Gamma\left(\nu + 3\right)} - \frac{x^{2}}{4 \Gamma\left(\nu + 2\right)} + \frac{1}{\Gamma\left(\nu + 1\right)}\right)}{\sin{\left(\pi \left(\nu + 1\right) \right)}} + \frac{\left(\frac{x}{2}\right)^{\nu + 1} \left(\left(\frac{x}{2}\right)^{\nu} \left(- \frac{x^{6}}{384 \Gamma\left(\nu + 4\right)} + \frac{x^{4}}{32 \Gamma\left(\nu + 3\right)} - \frac{x^{2}}{4 \Gamma\left(\nu + 2\right)} + \frac{1}{\Gamma\left(\nu + 1\right)}\right) \cos{\left(\pi \nu \right)} - \left(\frac{x}{2}\right)^{- \nu} \left(- \frac{x^{6}}{384 \Gamma\left(4 - \nu\right)} + \frac{x^{4}}{32 \Gamma\left(3 - \nu\right)} - \frac{x^{2}}{4 \Gamma\left(2 - \nu\right)} + \frac{1}{\Gamma\left(1 - \nu\right)}\right)\right) \left(- \frac{x^{6}}{384 \Gamma\left(\nu + 5\right)} + \frac{x^{4}}{32 \Gamma\left(\nu + 4\right)} - \frac{x^{2}}{4 \Gamma\left(\nu + 3\right)} + \frac{1}{\Gamma\left(\nu + 2\right)}\right)}{\sin{\left(\pi \nu \right)}}$
In [29]:
W3 = W2.expand(trig=True) + O(x**ordnung)
W3
Out[29]:
$\displaystyle - \frac{2}{x \sin{\left(\pi \nu \right)} \Gamma\left(- \nu\right) \Gamma\left(\nu + 1\right)} + \frac{x}{2 \sin{\left(\pi \nu \right)} \Gamma\left(- \nu\right) \Gamma\left(\nu + 2\right)} + \frac{x}{2 \sin{\left(\pi \nu \right)} \Gamma\left(1 - \nu\right) \Gamma\left(\nu + 1\right)} - \frac{x}{2 \sin{\left(\pi \nu \right)} \Gamma\left(1 - \nu\right) \Gamma\left(\nu + 2\right)} - \frac{x^{3}}{16 \sin{\left(\pi \nu \right)} \Gamma\left(- \nu\right) \Gamma\left(\nu + 3\right)} - \frac{x^{3}}{8 \sin{\left(\pi \nu \right)} \Gamma\left(1 - \nu\right) \Gamma\left(\nu + 2\right)} + \frac{x^{3}}{8 \sin{\left(\pi \nu \right)} \Gamma\left(1 - \nu\right) \Gamma\left(\nu + 3\right)} - \frac{x^{3}}{16 \sin{\left(\pi \nu \right)} \Gamma\left(2 - \nu\right) \Gamma\left(\nu + 1\right)} + \frac{x^{3}}{8 \sin{\left(\pi \nu \right)} \Gamma\left(2 - \nu\right) \Gamma\left(\nu + 2\right)} + O\left(x^{4}\right)$
In [30]:
simplify(W3.removeO())
Out[30]:
$\displaystyle \frac{2}{\pi x}$

Im richtigen Leben müsste man z.B. mit ordnung=8 und ordnung=20 wiederholen

Gekoppelte Pendel¶

In [31]:
y = Function('y')
w = Function('w')
t = Symbol('t', positive=True)
In [32]:
a = Rational(1,7)
b = Rational(1,6)
In [33]:
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
Out[33]:
$\displaystyle \left\{\frac{d^{2}}{d t^{2}} w{\left(t \right)} = - \frac{7 w{\left(t \right)}}{6} + \frac{y{\left(t \right)}}{6}, \frac{d^{2}}{d t^{2}} y{\left(t \right)} = \frac{w{\left(t \right)}}{7} - \frac{8 y{\left(t \right)}}{7}\right\}$
In [34]:
lsg = dsolve(dgs)
lsg
Out[34]:
$\displaystyle \left[ w{\left(t \right)} = - \frac{7 \sqrt{2310} C_{1} \sin{\left(\frac{\sqrt{2310} t}{42} \right)}}{330} - \frac{7 \sqrt{2310} C_{2} \cos{\left(\frac{\sqrt{2310} t}{42} \right)}}{330} + C_{3} \sin{\left(t \right)} + C_{4} \cos{\left(t \right)}, \ y{\left(t \right)} = \frac{\sqrt{2310} C_{1} \sin{\left(\frac{\sqrt{2310} t}{42} \right)}}{55} + \frac{\sqrt{2310} C_{2} \cos{\left(\frac{\sqrt{2310} t}{42} \right)}}{55} + C_{3} \sin{\left(t \right)} + C_{4} \cos{\left(t \right)}\right]$
In [35]:
phi_c = lsg[0].rhs
psi_c = lsg[1].rhs

Bestimme reelles Fundamentalsystem

In [36]:
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)
In [37]:
ersetzungen
Out[37]:
$\displaystyle \left[ \left\{ C_{1} : 1, \ C_{2} : 0, \ C_{3} : 0, \ C_{4} : 0\right\}, \ \left\{ C_{1} : 0, \ C_{2} : 1, \ C_{3} : 0, \ C_{4} : 0\right\}, \ \left\{ C_{1} : 0, \ C_{2} : 0, \ C_{3} : 1, \ C_{4} : 0\right\}, \ \left\{ C_{1} : 0, \ C_{2} : 0, \ C_{3} : 0, \ C_{4} : 1\right\}\right]$
In [38]:
fus = []
for s in ersetzungen:
    fus.append([phi_c.subs(s), 
                psi_c.subs(s)])
In [39]:
fus
Out[39]:
$\displaystyle \left[ \left[ - \frac{7 \sqrt{2310} \sin{\left(\frac{\sqrt{2310} t}{42} \right)}}{330}, \ \frac{\sqrt{2310} \sin{\left(\frac{\sqrt{2310} t}{42} \right)}}{55}\right], \ \left[ - \frac{7 \sqrt{2310} \cos{\left(\frac{\sqrt{2310} t}{42} \right)}}{330}, \ \frac{\sqrt{2310} \cos{\left(\frac{\sqrt{2310} t}{42} \right)}}{55}\right], \ \left[ \sin{\left(t \right)}, \ \sin{\left(t \right)}\right], \ \left[ \cos{\left(t \right)}, \ \cos{\left(t \right)}\right]\right]$
In [40]:
f = fus[0][0] + fus[2][0]
g = fus[0][1] + fus[2][1]
f, g
Out[40]:
$\displaystyle \left( \sin{\left(t \right)} - \frac{7 \sqrt{2310} \sin{\left(\frac{\sqrt{2310} t}{42} \right)}}{330}, \ \sin{\left(t \right)} + \frac{\sqrt{2310} \sin{\left(\frac{\sqrt{2310} t}{42} \right)}}{55}\right)$
In [41]:
plot(f,g+3, (t,-100,100));