Lektion 7



In [1]:

    
from sympy import *
init_printing()
import numpy as np
import matplotlib.pyplot as plt

interaktiv



In [ ]:

    
%matplotlib notebook

für den Druck



In [2]:

    
%matplotlib inline

Integration



In [3]:

    
x = Symbol('x')



In [4]:

    
f = sin(x)/x



In [5]:

    
I1 = Integral(f, x)
I1









    Out[5]:





$$\int \frac{1}{x} \sin{\left (x \right )}\, dx$$



In [7]:

    
I1.doit()









    Out[7]:





$$\operatorname{Si}{\left (x \right )}$$



In [8]:

    
plot(Si(x), (x, -8, 8));



In [9]:

    
Si(6.)









    Out[9]:





$$1.42468755128051$$



In [10]:

    
Integral(sin(x)/x, (x,0,6)).n()









    Out[10]:





$$1.42468755128051$$



In [11]:

    
I3 = Integral(exp(-x**2), x)
I3









    Out[11]:





$$\int e^{- x^{2}}\, dx$$



In [12]:

    
I3.doit()









    Out[12]:





$$\frac{\sqrt{\pi}}{2} \operatorname{erf}{\left (x \right )}$$



In [13]:

    
plot(erf(x), (x, -3, 3));



In [14]:

    
f = (1-x**2)**Rational(-3,2)



In [15]:

    
I4 = Integral(f, x)
I4









    Out[15]:





$$\int \frac{1}{\left(- x^{2} + 1\right)^{\frac{3}{2}}}\, dx$$



In [16]:

    
tmp = I4.doit()
tmp









    Out[16]:





$$\begin{cases} \frac{1}{\sqrt{-1 + \frac{1}{x^{2}}}} & \text{for}\: \left|{\frac{1}{x^{2}}}\right| > 1 \\- \frac{i}{\sqrt{1 - \frac{1}{x^{2}}}} & \text{otherwise} \end{cases}$$



In [17]:

    
F = tmp.args[0][0]
F









    Out[17]:





$$\frac{1}{\sqrt{-1 + \frac{1}{x^{2}}}}$$



In [18]:

    
F.diff(x) == f









    Out[18]:





False



In [19]:

    
F.diff(x)









    Out[19]:





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



In [23]:

    
(F.diff(x) - f).radsimp().expand().ratsimp()









    Out[23]:





$$\frac{x \sqrt{-1 + \frac{1}{x^{2}}} - \sqrt{- x^{2} + 1}}{x^{4} - 2 x^{2} + 1}$$



In [24]:

    
xn = np.linspace(-.8, .8)
fn = lambdify(x, f)
Fn = lambdify(x, F, 'numpy')
dFn = lambdify(x, F.diff(x), 'numpy')
plt.plot(xn, fn(xn), label='f')
plt.plot(xn, Fn(xn), label='F')
plt.plot(xn, dFn(xn), label='$\\frac{\\partial F}{\\partial x}$')
plt.legend(loc='lower right');

Die Stammfunktion ist falsch. Wir waren aber gewarnt.



In [25]:

    
I5 = Integral(f, (x, Rational(1,5), Rational(1,2)))
I5









    Out[25]:





$$\int_{\frac{1}{5}}^{\frac{1}{2}} \frac{1}{\left(- x^{2} + 1\right)^{\frac{3}{2}}}\, dx$$



In [26]:

    
I5.doit()  # AttribureError









    



---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-26-40452399f355> in <module>()
----> 1 I5.doit()  # AttribureError

C:\Users\braun\Miniconda3\lib\site-packages\sympy\integrals\integrals.py in doit(self, **hints)
    551                                        for f in integrals])
    552                         try:
--> 553                             evalued = Add(*others)._eval_interval(x, a, b)
    554                             function = uneval + evalued
    555                         except NotImplementedError:

C:\Users\braun\Miniconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in _eval_interval(self, sym, a, b)
    270 
    271         # Determine what intervals the expr,cond pairs affect.
--> 272         int_expr = self._sort_expr_cond(sym, a, b)
    273 
    274         # Finally run through the intervals and sum the evaluation.

C:\Users\braun\Miniconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in _sort_expr_cond(self, sym, a, b, targetcond)
    341                 if sym not in [cond.lts, cond.gts]:
    342                     cond = _solve_inequality(cond, sym)
--> 343                 lower, upper = cond.lts, cond.gts  # part 1: initialize with givens
    344                 if cond.lts == sym:                # part 1a: expand the side ...
    345                     lower = S.NegativeInfinity   # e.g. x <= 0 ---> -oo <= 0

AttributeError: 'Or' object has no attribute 'lts'



In [27]:

    
y = Symbol('y', real=True)



In [28]:

    
g = sin(x) * exp(-x*y)
I6 = Integral(g, x)
I6









    Out[28]:





$$\int e^{- x y} \sin{\left (x \right )}\, dx$$



In [29]:

    
G = I6.doit()
G









    Out[29]:





$$- \frac{y \sin{\left (x \right )}}{y^{2} e^{x y} + e^{x y}} - \frac{\cos{\left (x \right )}}{y^{2} e^{x y} + e^{x y}}$$



In [31]:

    
G.diff(x)# == g









    Out[31]:





$$- \frac{y \cos{\left (x \right )}}{y^{2} e^{x y} + e^{x y}} - \frac{y \sin{\left (x \right )}}{\left(y^{2} e^{x y} + e^{x y}\right)^{2}} \left(- y^{3} e^{x y} - y e^{x y}\right) + \frac{\sin{\left (x \right )}}{y^{2} e^{x y} + e^{x y}} - \frac{\left(- y^{3} e^{x y} - y e^{x y}\right) \cos{\left (x \right )}}{\left(y^{2} e^{x y} + e^{x y}\right)^{2}}$$



In [32]:

    
G.diff(x).simplify() == g









    Out[32]:





True

Das Taylorpolynom



In [33]:

    
x = Symbol('x')
f = sqrt(1+x)/sqrt(1-x**2)
f









    Out[33]:





$$\frac{\sqrt{x + 1}}{\sqrt{- x^{2} + 1}}$$



In [34]:

    
t = f.series(x, 0, 8)
t









    Out[34]:





$$1 + \frac{x}{2} + \frac{3 x^{2}}{8} + \frac{5 x^{3}}{16} + \frac{35 x^{4}}{128} + \frac{63 x^{5}}{256} + \frac{231 x^{6}}{1024} + \frac{429 x^{7}}{2048} + \mathcal{O}\left(x^{8}\right)$$



In [35]:

    
xn = np.linspace(-.8, .8)
fn = lambdify(x, f, 'numpy')
plt.plot(xn, fn(xn), label='f')
for n in range(2, 5):   # n=1 führt auf ValueError
    tt = t + O(x**n, (x,0))
    tn = lambdify(x, tt.removeO(), 'numpy')
    plt.plot(xn, tn(xn), label='n='+str(n))
plt.legend(loc='upper left');

Dasselbe für den Arcustangens



In [36]:

    
xn = np.linspace(-2, 2, 200)
plt.plot(xn, np.arctan(xn), label='arctan')
for n in [4, 20, 60]:
    t = atan(x).series(x, 0, n)
    tn = lambdify(x, t.removeO(), 'numpy')
    plt.plot(xn, tn(xn), label='n='+str(n))
plt.axis(ymin=-2, ymax=2)
plt.legend(loc='upper center');

Allgemeinere Reihenentwicklungen



In [45]:

    
t = {}
t[1] = atan(x).series(x, -oo, 12)
t[1]









    Out[45]:





$$\frac{1}{11 x^{11}} - \frac{1}{9 x^{9}} + \frac{1}{7 x^{7}} - \frac{1}{5 x^{5}} + \frac{1}{3 x^{3}} - \frac{1}{x} - \frac{\pi}{2} + \mathcal{O}\left(\frac{1}{x^{12}}; x\rightarrow\infty\right)$$



In [47]:

    
t[2] = atan(x).series(x, 0, 12)
t[2]









    Out[47]:





$$x - \frac{x^{3}}{3} + \frac{x^{5}}{5} - \frac{x^{7}}{7} + \frac{x^{9}}{9} - \frac{x^{11}}{11} + \mathcal{O}\left(x^{12}\right)$$



In [48]:

    
t[3] = atan(x).series(x, oo, 12)
t[3]









    Out[48]:





$$\frac{1}{11 x^{11}} - \frac{1}{9 x^{9}} + \frac{1}{7 x^{7}} - \frac{1}{5 x^{5}} + \frac{1}{3 x^{3}} - \frac{1}{x} + \frac{\pi}{2} + \mathcal{O}\left(\frac{1}{x^{12}}; x\rightarrow\infty\right)$$



In [49]:

    
xn = np.linspace(-3, 3, 200)
for nr, tt in t.items():
    tn = lambdify(x, tt.removeO())
    plt.plot(xn, tn(xn))
plt.axis(ymin=-2, ymax=2);

Grenzwertbestimmung mittels Reihenentwicklung



In [50]:

    
x = Symbol('x')



In [51]:

    
f = 1 - cos(x**2)
g = x*(x-sin(x))
L = Limit(f/g, x, 0)
L









    Out[51]:





$$\lim_{x \to 0^+}\left(\frac{- \cos{\left (x^{2} \right )} + 1}{x \left(x - \sin{\left (x \right )}\right)}\right)$$



In [52]:

    
L.doit()









    Out[52]:





$$3$$



In [53]:

    
fr = f.series(x, 0, 10)
fr









    Out[53]:





$$\frac{x^{4}}{2} - \frac{x^{8}}{24} + \mathcal{O}\left(x^{10}\right)$$



In [54]:

    
gr = g.series(x, 0, 10)
gr









    Out[54]:





$$\frac{x^{4}}{6} - \frac{x^{6}}{120} + \frac{x^{8}}{5040} + \mathcal{O}\left(x^{10}\right)$$



In [55]:

    
t1 = fr.removeO() / gr.removeO()
t1









    Out[55]:





$$\frac{- \frac{x^{8}}{24} + \frac{x^{4}}{2}}{\frac{x^{8}}{5040} - \frac{x^{6}}{120} + \frac{x^{4}}{6}}$$



In [56]:

    
n, z = numer(t1), denom(t1)
t2 = (n/x**4).expand() / (z/x**4).expand()
t2









    Out[56]:





$$\frac{- \frac{x^{4}}{24} + \frac{1}{2}}{\frac{x^{4}}{5040} - \frac{x^{2}}{120} + \frac{1}{6}}$$



In [57]:

    
t2.subs(x, 0)









    Out[57]:





$$3$$

Vektoren und Matrizen



In [58]:

    
x = Matrix([1,2,3])
x









    Out[58]:





$$\left[\begin{matrix}1\\2\\3\end{matrix}\right]$$



In [59]:

    
y = Matrix(1,3,[4,5,6])
y









    Out[59]:





$$\left[\begin{matrix}4 & 5 & 6\end{matrix}\right]$$



In [60]:

    
y*x









    Out[60]:





$$\left[\begin{matrix}32\end{matrix}\right]$$



In [61]:

    
x*y









    Out[61]:





$$\left[\begin{matrix}4 & 5 & 6\\8 & 10 & 12\\12 & 15 & 18\end{matrix}\right]$$



In [62]:

    
A = Matrix(3,3,range(1,10))
A









    Out[62]:





$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$$



In [63]:

    
A[1,1]









    Out[63]:





$$5$$



In [64]:

    
A*x









    Out[64]:





$$\left[\begin{matrix}14\\32\\50\end{matrix}\right]$$



In [65]:

    
A*y  # ShapeError









    



---------------------------------------------------------------------------
ShapeError                                Traceback (most recent call last)
<ipython-input-65-b68ff40d08ab> in <module>()
----> 1 A*y  # ShapeError

C:\Users\braun\Miniconda3\lib\site-packages\sympy\core\decorators.py in binary_op_wrapper(self, other)
    116                     else:
    117                         return f(self)
--> 118             return func(self, other)
    119         return binary_op_wrapper
    120     return priority_decorator

C:\Users\braun\Miniconda3\lib\site-packages\sympy\matrices\dense.py in __mul__(self, other)
    553     @call_highest_priority('__rmul__')
    554     def __mul__(self, other):
--> 555         return super(DenseMatrix, self).__mul__(_force_mutable(other))
    556 
    557     @call_highest_priority('__rmul__')

C:\Users\braun\Miniconda3\lib\site-packages\sympy\matrices\matrices.py in __mul__(self, other)
    502             B = other
    503             if A.cols != B.rows:
--> 504                 raise ShapeError("Matrices size mismatch.")
    505             if A.cols == 0:
    506                 return classof(A, B)._new(A.rows, B.cols, lambda i, j: 0)

ShapeError: Matrices size mismatch.



In [66]:

    
y*A









    Out[66]:





$$\left[\begin{matrix}66 & 81 & 96\end{matrix}\right]$$



In [67]:

    
B = Matrix(3,3,[1,0,1,-1,-1,-1,0,1,0])
B









    Out[67]:





$$\left[\begin{matrix}1 & 0 & 1\\-1 & -1 & -1\\0 & 1 & 0\end{matrix}\right]$$



In [68]:

    
A * B









    Out[68]:





$$\left[\begin{matrix}-1 & 1 & -1\\-1 & 1 & -1\\-1 & 1 & -1\end{matrix}\right]$$



In [70]:

    
np.array(A) * np.array(B)









    Out[70]:





array([[1, 0, 3],
       [-4, -5, -6],
       [0, 8, 0]], dtype=object)



In [71]:

    
A*B -B*A









    Out[71]:





$$\left[\begin{matrix}-9 & -9 & -13\\11 & 16 & 17\\-5 & -4 & -7\end{matrix}\right]$$



In [72]:

    
eye(5)









    Out[72]:





$$\left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 1 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 1\end{matrix}\right]$$



In [74]:

    
C = A + eye(3)



In [75]:

    
A.det()









    Out[75]:





$$0$$



In [76]:

    
C.det()









    Out[76]:





$$-2$$



In [77]:

    
C**(-1)









    Out[77]:





$$\left[\begin{matrix}-6 & -2 & 3\\-1 & \frac{1}{2} & 0\\5 & 1 & -2\end{matrix}\right]$$



In [78]:

    
C*C**(-1)









    Out[78]:





$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$



In [79]:

    
C**(-1)*C









    Out[79]:





$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$



In [80]:

    
A.T  # Transponierte









    Out[80]:





$$\left[\begin{matrix}1 & 4 & 7\\2 & 5 & 8\\3 & 6 & 9\end{matrix}\right]$$

Manipulation von Matrizen



In [81]:

    
Matrix([A,B])









    Out[81]:





$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\\1 & 0 & 1\\-1 & -1 & -1\\0 & 1 & 0\end{matrix}\right]$$



In [82]:

    
Matrix([A.T, B.T]).T









    Out[82]:





$$\left[\begin{matrix}1 & 2 & 3 & 1 & 0 & 1\\4 & 5 & 6 & -1 & -1 & -1\\7 & 8 & 9 & 0 & 1 & 0\end{matrix}\right]$$



In [85]:

    
A.reshape(1, len(A))#.reshape(3,3)   # flatten a matrix









    Out[85]:





$$\left[\begin{matrix}1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9\end{matrix}\right]$$



In [86]:

    
A









    Out[86]:





$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$$



In [87]:

    
A[0:2, 1:3]









    Out[87]:





$$\left[\begin{matrix}2 & 3\\5 & 6\end{matrix}\right]$$



In [88]:

    
A[[0,2], [0,2]]









    Out[88]:





$$\left[\begin{matrix}1 & 3\\7 & 9\end{matrix}\right]$$

Matrix aus Bildungsvorschrift



In [89]:

    
def element(i,j):
    return 1 + 3*i + j



In [90]:

    
E = Matrix(3,3, element)
E









    Out[90]:





$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\end{matrix}\right]$$



In [91]:

    
def hilbert_element(i,j):
    return 1/(i+j+1)



In [92]:

    
hilbert = Matrix(5,5,hilbert_element)
hilbert









    Out[92]:





$$\left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7}\\\frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8}\\\frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9}\end{matrix}\right]$$



In [93]:

    
hilbert.det()









    Out[93]:





$$\frac{1}{266716800000}$$



In [94]:

    
hilbert**(-1)









    Out[94]:





$$\left[\begin{matrix}25 & -300 & 1050 & -1400 & 630\\-300 & 4800 & -18900 & 26880 & -12600\\1050 & -18900 & 79380 & -117600 & 56700\\-1400 & 26880 & -117600 & 179200 & -88200\\630 & -12600 & 56700 & -88200 & 44100\end{matrix}\right]$$



In [95]:

    
hilbert = Matrix(60, 60, hilbert_element)



In [96]:

    
inv_hilbert = hilbert**(-1)



In [97]:

    
inv_hilbert[29,40]









    Out[97]:





$$-13106188877111499639796065719840267765492147640133102288707468863670920137401990144000$$



In [98]:

    
def f_hilbert_element(i,j):
    return float(hilbert_element(i,j))

f_hilbert = Matrix(60, 60, f_hilbert_element)



In [99]:

    
inv_f_hilbert = f_hilbert**(-1)



In [100]:

    
inv_f_hilbert[29,40]









    Out[100]:





$$-1.31643790355102 \cdot 10^{16}$$



In [101]:

    
float(inv_hilbert[29,40])









    Out[101]:





$$-1.31061888771115e+85$$



In [102]:

    
float(hilbert.det())









    Out[102]:





$$0.0$$



In [103]:

    
hd = hilbert.det()
hd









    Out[103]:





$$\frac{1}{144595578563805514086158249720372269159002085903382430132201033126808981249036023606637008387023829599333666316075604948312115145411863969712946077100121030682846678058070583611393911862619093902501386471375819221747631611277681985408821087315622109090337413025602702195139847390738623762934604821752342674239542453237619108127011300097716120691693898949607679638036065375186860306076222363870724858462957501450975060651452580879958066769635565836268727971708263260406956936283896739981777305823220951557849587117301546890968416334028835991158624545856654613088266798404474173809571065306900215468090502720547807378701429903835996074831301035858455007298927233971501116116274485491900273784036359713781894831240582761862528401276169429713774952186680229710780217577488661819819689511945891416643911424747998815748756279920991164032891287911502462152209429901593587296320296206234744776753094385143872208007096501328349984197696099481287704088458246345934951881581833852029911538394663184463632483139054401713155038298109933169825010153712197777995718512077459470886069933743445187868699045328155791603974619737440747367977499216895232157255157009596567132433587296576330286987646097003836387218585334927870388584495533920436544777246123572137496212023370398668473560704584359848492007996040151478841677215515073951567362180563147602372933568571677649505118065555493148990970392349507555863857592128666495835904521226456219958978767742318208524954164146330757260249648322351081668227716769899597203478929347736483260126579844187151269219354808104825593392288847315732218934859633534872295402195352956429504653738061031861366606590885485118028082106401430020143303222325054385500356570233107440548265674758703842629685083189043850626308746522354010231269905840690203495093487771272862040284657239551505054851451998499232307307321353832093108718887300229729038702644191598020531876236248898677746718045468357283324738175869625155835037278573318947974079408327907156769659049187005858872226622838121073847595276697600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000}$$



In [104]:

    
mpmath.mpf(hd)









    Out[104]:





mpf('6.9158407880274928e-2121')



In [ ]: