Lektion 13



In [3]:

    
from sympy import *
init_printing()
from IPython.display import display

Besselfunktionen



In [4]:

    
x = Symbol('x')
y = Function('y')
dgl = Eq(y(x).diff(x, 2), -1/x*y(x).diff(x) + 1/x**2*y(x) +4*y(x))
dgl









    Out[4]:





$$\frac{d^{2}}{d x^{2}}  y{\left (x \right )} = 4 y{\left (x \right )} - \frac{1}{x} \frac{d}{d x} y{\left (x \right )} + \frac{1}{x^{2}} y{\left (x \right )}$$



In [6]:

    
#dsolve(dgl)   # NotImplementedError



In [28]:

    
#N = 8
N=18



In [29]:

    
a = [Symbol('a'+str(j)) for j in range(N)]



In [30]:

    
n = Symbol('n')



In [31]:

    
ys = sum([a[j]*x**j for j in range(N)]) 
ys









    Out[31]:





$$a_{0} + a_{1} x + a_{10} x^{10} + a_{11} x^{11} + a_{12} x^{12} + a_{13} x^{13} + a_{14} x^{14} + a_{15} x^{15} + a_{16} x^{16} + a_{17} x^{17} + a_{2} x^{2} + a_{3} x^{3} + a_{4} x^{4} + a_{5} x^{5} + a_{6} x^{6} + a_{7} x^{7} + a_{8} x^{8} + a_{9} x^{9}$$



In [32]:

    
gl = dgl.subs(y(x), ys).doit()
gl









    Out[32]:





$$2 \left(45 a_{10} x^{8} + 55 a_{11} x^{9} + 66 a_{12} x^{10} + 78 a_{13} x^{11} + 91 a_{14} x^{12} + 105 a_{15} x^{13} + 120 a_{16} x^{14} + 136 a_{17} x^{15} + a_{2} + 3 a_{3} x + 6 a_{4} x^{2} + 10 a_{5} x^{3} + 15 a_{6} x^{4} + 21 a_{7} x^{5} + 28 a_{8} x^{6} + 36 a_{9} x^{7}\right) = 4 a_{0} + 4 a_{1} x + 4 a_{10} x^{10} + 4 a_{11} x^{11} + 4 a_{12} x^{12} + 4 a_{13} x^{13} + 4 a_{14} x^{14} + 4 a_{15} x^{15} + 4 a_{16} x^{16} + 4 a_{17} x^{17} + 4 a_{2} x^{2} + 4 a_{3} x^{3} + 4 a_{4} x^{4} + 4 a_{5} x^{5} + 4 a_{6} x^{6} + 4 a_{7} x^{7} + 4 a_{8} x^{8} + 4 a_{9} x^{9} - \frac{1}{x} \left(a_{1} + 10 a_{10} x^{9} + 11 a_{11} x^{10} + 12 a_{12} x^{11} + 13 a_{13} x^{12} + 14 a_{14} x^{13} + 15 a_{15} x^{14} + 16 a_{16} x^{15} + 17 a_{17} x^{16} + 2 a_{2} x + 3 a_{3} x^{2} + 4 a_{4} x^{3} + 5 a_{5} x^{4} + 6 a_{6} x^{5} + 7 a_{7} x^{6} + 8 a_{8} x^{7} + 9 a_{9} x^{8}\right) + \frac{1}{x^{2}} \left(a_{0} + a_{1} x + a_{10} x^{10} + a_{11} x^{11} + a_{12} x^{12} + a_{13} x^{13} + a_{14} x^{14} + a_{15} x^{15} + a_{16} x^{16} + a_{17} x^{17} + a_{2} x^{2} + a_{3} x^{3} + a_{4} x^{4} + a_{5} x^{5} + a_{6} x^{6} + a_{7} x^{7} + a_{8} x^{8} + a_{9} x^{9}\right)$$



In [33]:

    
p1 = (gl.lhs - gl.rhs).expand()
p1









    Out[33]:





$$- 4 a_{0} - \frac{a_{0}}{x^{2}} - 4 a_{1} x - 4 a_{10} x^{10} + 99 a_{10} x^{8} - 4 a_{11} x^{11} + 120 a_{11} x^{9} - 4 a_{12} x^{12} + 143 a_{12} x^{10} - 4 a_{13} x^{13} + 168 a_{13} x^{11} - 4 a_{14} x^{14} + 195 a_{14} x^{12} - 4 a_{15} x^{15} + 224 a_{15} x^{13} - 4 a_{16} x^{16} + 255 a_{16} x^{14} - 4 a_{17} x^{17} + 288 a_{17} x^{15} - 4 a_{2} x^{2} + 3 a_{2} - 4 a_{3} x^{3} + 8 a_{3} x - 4 a_{4} x^{4} + 15 a_{4} x^{2} - 4 a_{5} x^{5} + 24 a_{5} x^{3} - 4 a_{6} x^{6} + 35 a_{6} x^{4} - 4 a_{7} x^{7} + 48 a_{7} x^{5} - 4 a_{8} x^{8} + 63 a_{8} x^{6} - 4 a_{9} x^{9} + 80 a_{9} x^{7}$$



In [34]:

    
p1.coeff(x**(-2))









    Out[34]:





$$- a_{0}$$



In [35]:

    
p1.coeff(x**(-1))









    Out[35]:





$$0$$



In [36]:

    
p1.coeff(x, 1)









    Out[36]:





$$- 4 a_{1} + 8 a_{3}$$



In [37]:

    
gls = []
for j in range(N+1):
    glg = Eq(p1.coeff(x, j-2), 0)
    if glg != True:
        gls.append(glg)
gls









    Out[37]:





$$\left [ - a_{0} = 0, \quad - 4 a_{0} + 3 a_{2} = 0, \quad - 4 a_{1} + 8 a_{3} = 0, \quad - 4 a_{2} + 15 a_{4} = 0, \quad - 4 a_{3} + 24 a_{5} = 0, \quad - 4 a_{4} + 35 a_{6} = 0, \quad - 4 a_{5} + 48 a_{7} = 0, \quad - 4 a_{6} + 63 a_{8} = 0, \quad - 4 a_{7} + 80 a_{9} = 0, \quad 99 a_{10} - 4 a_{8} = 0, \quad 120 a_{11} - 4 a_{9} = 0, \quad - 4 a_{10} + 143 a_{12} = 0, \quad - 4 a_{11} + 168 a_{13} = 0, \quad - 4 a_{12} + 195 a_{14} = 0, \quad - 4 a_{13} + 224 a_{15} = 0, \quad - 4 a_{14} + 255 a_{16} = 0, \quad - 4 a_{15} + 288 a_{17} = 0, \quad - 4 a_{16} = 0\right ]$$



In [23]:

    
#solve(gls) #NotImplementedError









    Out[23]:





$$\left \{ a_{0} : 0, \quad a_{1} : 144 a_{7}, \quad a_{2} : 0, \quad a_{3} : 72 a_{7}, \quad a_{4} : 0, \quad a_{5} : 12 a_{7}, \quad a_{6} : 0\right \}$$



In [38]:

    
Lsg = solve(gls[:-1])
Lsg









    Out[38]:





$$\left \{ a_{0} : 0, \quad a_{1} : 2880 a_{9}, \quad a_{10} : 0, \quad a_{11} : \frac{a_{9}}{30}, \quad a_{12} : 0, \quad a_{13} : \frac{a_{9}}{1260}, \quad a_{14} : 0, \quad a_{15} : \frac{a_{9}}{70560}, \quad a_{16} : 0, \quad a_{17} : \frac{a_{9}}{5080320}, \quad a_{2} : 0, \quad a_{3} : 1440 a_{9}, \quad a_{4} : 0, \quad a_{5} : 240 a_{9}, \quad a_{6} : 0, \quad a_{7} : 20 a_{9}, \quad a_{8} : 0\right \}$$

Die ungeraden a werden rückwärts gelöst. Das ist verwirrend.



In [39]:

    
var = a.copy()   # böse Falle
del var[1]
var









    Out[39]:





$$\left [ a_{0}, \quad a_{2}, \quad a_{3}, \quad a_{4}, \quad a_{5}, \quad a_{6}, \quad a_{7}, \quad a_{8}, \quad a_{9}, \quad a_{10}, \quad a_{11}, \quad a_{12}, \quad a_{13}, \quad a_{14}, \quad a_{15}, \quad a_{16}, \quad a_{17}\right ]$$



In [40]:

    
Lsg = solve(gls[:-1], var)
Lsg









    Out[40]:





$$\left \{ a_{0} : 0, \quad a_{10} : 0, \quad a_{11} : \frac{a_{1}}{86400}, \quad a_{12} : 0, \quad a_{13} : \frac{a_{1}}{3628800}, \quad a_{14} : 0, \quad a_{15} : \frac{a_{1}}{203212800}, \quad a_{16} : 0, \quad a_{17} : \frac{a_{1}}{14631321600}, \quad a_{2} : 0, \quad a_{3} : \frac{a_{1}}{2}, \quad a_{4} : 0, \quad a_{5} : \frac{a_{1}}{12}, \quad a_{6} : 0, \quad a_{7} : \frac{a_{1}}{144}, \quad a_{8} : 0, \quad a_{9} : \frac{a_{1}}{2880}\right \}$$

Wir hatten das beim ersten Mal mit $N=8$ gemacht. Das sind zu wenige Daten. Jetzt noch Mal mit $N=18$.

Aus Gründen, die ich nicht verstehe, muss man den Kernel zurücksetzen, bevor man mit dem neuen $N$ startet.



In [42]:

    
#raise Unterbrechung



In [43]:

    
Lsg[a[1]] = a[1]



In [44]:

    
q = [Lsg[a[2*j+3]]/Lsg[a[2*j+1]] for j in range(int(N/2)-2)]
display(q)









    





$$\left [ \frac{1}{2}, \quad \frac{1}{6}, \quad \frac{1}{12}, \quad \frac{1}{20}, \quad \frac{1}{30}, \quad \frac{1}{42}, \quad \frac{1}{56}\right ]$$



In [45]:

    
liste = []
for j in range(int(N/2-2)):
    m = Lsg[a[2*j+1]]/Lsg[a[2*j+3]]
    liste.append(m/(j+2))
liste









    Out[45]:





$$\left [ 1, \quad 2, \quad 3, \quad 4, \quad 5, \quad 6, \quad 7\right ]$$

Also $$ \frac{a_{2j+3}}{a_{2j+1}} = \frac1{(j+1)(j+2)} $$

Das bedeutet $$ a_{2n+3} = \prod_{j=0}^n \frac1{(j+1)(j+2)} a_1 = \frac{a_1}{(n+1)!(n+2)!}. $$

Probe



In [46]:

    
for j in range(int(N/2)-2):
    display((Lsg[a[2*j+1]], a[1]/factorial(j)/factorial(j+1)))









    





$$\left ( a_{1}, \quad a_{1}\right )$$






    





$$\left ( \frac{a_{1}}{2}, \quad \frac{a_{1}}{2}\right )$$






    





$$\left ( \frac{a_{1}}{12}, \quad \frac{a_{1}}{12}\right )$$






    





$$\left ( \frac{a_{1}}{144}, \quad \frac{a_{1}}{144}\right )$$






    





$$\left ( \frac{a_{1}}{2880}, \quad \frac{a_{1}}{2880}\right )$$






    





$$\left ( \frac{a_{1}}{86400}, \quad \frac{a_{1}}{86400}\right )$$






    





$$\left ( \frac{a_{1}}{3628800}, \quad \frac{a_{1}}{3628800}\right )$$



In [47]:

    
S1 = Sum(x**(2*n+1)/factorial(n)/factorial(n+1), (n,0,oo))
S1









    Out[47]:





$$\sum_{n=0}^{\infty} \frac{x^{2 n + 1}}{n! \left(n + 1\right)!}$$



In [48]:

    
u = S1.doit()
u









    Out[48]:





$$I_{1}\left(2 x\right)$$



In [49]:

    
srepr(u)









    Out[49]:





"besseli(Integer(1), Mul(Integer(2), Symbol('x')))"



In [51]:

    
besseli?

Eine zweite Lösung müsste man erhalten können, indem man einen Ansatz aus einer Potenzreihe und dem Produkt aus dem Logarithmus und einer Potenzreihe macht.

Das Reduktionsverfahren von d'Alembert führt auf ein schwieriges Integral.



In [52]:

    
tmp = dgl.subs(y(x), u).doit()
tmp









    Out[52]:





$$3 I_{1}\left(2 x\right) + I_{3}\left(2 x\right) = 4 I_{1}\left(2 x\right) - \frac{1}{x} \left(I_{0}\left(2 x\right) + I_{2}\left(2 x\right)\right) + \frac{I_{1}\left(2 x\right)}{x^{2}}$$

http://dlmf.nist.gov

$$ I_{\nu-1}(z) - I_{\nu+1}(z) = \frac{2\nu}z I_\nu(z) $$



In [53]:

    
(tmp.lhs - tmp.rhs).series(x, 0, 20)









    Out[53]:





$$\mathcal{O}\left(x^{20}\right)$$

Pattern matching



In [54]:

    
x = Symbol('x')



In [55]:

    
x1 = Wild('x1')



In [56]:

    
pattern = sin(x1)



In [57]:

    
a = sin(2*x+5)



In [58]:

    
m = a.match(pattern)
m









    Out[58]:





$$\left \{ x_{1} : 2 x + 5\right \}$$



In [59]:

    
b = 2*sin(x1/2)*cos(x1/2)
b









    Out[59]:





$$2 \sin{\left (\frac{x_{1}}{2} \right )} \cos{\left (\frac{x_{1}}{2} \right )}$$



In [60]:

    
b.subs(m)









    Out[60]:





$$2 \sin{\left (x + \frac{5}{2} \right )} \cos{\left (x + \frac{5}{2} \right )}$$



In [61]:

    
def expand_sin_x_halbe(term):
    x1 = Wild('x1')
    pattern = sin(x1)
    ersetzung = 2*sin(x1/2)*cos(x1/2)
    m = term.match(pattern)
    if m:
        return ersetzung.subs(m)
    else:
        return term



In [62]:

    
expand_sin_x_halbe(sin(x/2))









    Out[62]:





$$2 \sin{\left (\frac{x}{4} \right )} \cos{\left (\frac{x}{4} \right )}$$



In [63]:

    
series(expand_sin_x_halbe(sin(x)) - sin(x), x, 0, 20)









    Out[63]:





$$\mathcal{O}\left(x^{20}\right)$$



In [64]:

    
a = 2*sin(x)
expand_sin_x_halbe(a)









    Out[64]:





$$2 \sin{\left (x \right )}$$



In [65]:

    
a.is_Mul









    Out[65]:





True



In [66]:

    
a.args









    Out[66]:





$$\left ( 2, \quad \sin{\left (x \right )}\right )$$



In [69]:

    
def expand_sin_x_halbe(ausdr):
    ausdr = S(ausdr)
    x1 = Wild('x1')
    pattern = sin(x1)
    ersetzung = 2*sin(x1/2)*cos(x1/2)
    m = ausdr.match(pattern)
    if m:
        res = ersetzung.subs(m)
    elif ausdr.is_Mul:
        res = 1
        for term in ausdr.args:
            res = res * expand_sin_x_halbe(term)
    elif ausdr.is_Add:
        res = 0
        for term in ausdr.args:
            res = res + expand_sin_x_halbe(term)
    else:
        res = ausdr
    return res



In [70]:

    
expand_sin_x_halbe(sin(2*x))









    Out[70]:





$$2 \sin{\left (x \right )} \cos{\left (x \right )}$$



In [71]:

    
expand_sin_x_halbe(sin(2*x)/2)









    Out[71]:





$$\sin{\left (x \right )} \cos{\left (x \right )}$$



In [72]:

    
expand_sin_x_halbe(1+10*sin((x+1)**2))









    Out[72]:





$$20 \sin{\left (\frac{1}{2} \left(x + 1\right)^{2} \right )} \cos{\left (\frac{1}{2} \left(x + 1\right)^{2} \right )} + 1$$



In [73]:

    
ausdr = (1 + sin(x/2))**3



In [74]:

    
expand_sin_x_halbe(ausdr)









    Out[74]:





$$\left(\sin{\left (\frac{x}{2} \right )} + 1\right)^{3}$$



In [75]:

    
ausdr.is_Pow









    Out[75]:





True

Mit etwas Fleiß bekommt man ein vollständiges Ersetzungssystem



In [ ]: