Musterlösung: Aufgabe 1, Übungsblatt 7

Einführung in die Programmierung für Studierende der Physik

(Introduction to Programming for Physicists)

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main

(Sommersemester 2022)

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske

Frankfurt am Main 28.05.2022

Kurvendiskussion mittels eines Jupyter Notebooks in Verbindung mit einem C++

Aufgabe 1

Im Folgenden soll ein Teil einer Kurvendiskussion mittels eines eigenen Jupyter Notebooks in Verbindung mit einem C++ angefertigt werden. Gegeben sei die Funktion $f(x)$ $$ \begin{equation} f(x) = e^{-x^2/4} \cdot \hbox{sin}(x/10) \cdot \left( x^4 -q\,x^2 - 5 \right) \quad . \end{equation} $$

Nullstellensuche mit Jupyter und Sympy

Betrachten Sie sich die Funktion in ihrem gesamten Definitionsbereich ($x \in$ ℝ) und berechnen Sie für beliebige Parameter $q \in$ ℝ die Nullstellen der Funktion $f(x)=0$. Wie lauten die Nullstellenwerte für $q=9$ ?

Darstellung von $f(x)$, $f^\prime(x)$ und $f^{\prime\prime}(x)$ mittels Jupyter

Stellen Sie nun die Funktion $f(x)$ und die erste und zweite Ableitung der Funktion im Teilintervall $x \in [-7,7]$ grafisch in drei x-y-Diagrammen dar (benutzen Sie dafür wieder den festen Parameterwert $q=9$). Versuchen Sie die Maxima der Funktion (mit $q=9$) mittels Sympy zu bestimmen.

C++ Programm zur Bestimmung der Hochpunkte der Funktion

Berechnen Sie die x-Werte der Hochpunkte der Funktion $f(x)$ mittels der Methode der Bisektion (Intervallhalbierungsverfahren) unter Verwendung eines C++ Programms.

Nullstellensuche mit Jupyter und Sympy

Wir definieren die Funktion und die Gleichung für die Nullstellensuche lautet:

In [1]:
from sympy import *
init_printing()
In [2]:
x = symbols('x')
q = symbols('q')
f = exp(-x**2/4)*sin(x/10)*(x**4 -q*x**2 - 5)
Gl_Nullstelle = Eq(f, 0)
Gl_Nullstelle
Out[2]:
$$\left(- q x^{2} + x^{4} - 5\right) e^{- \frac{x^{2}}{4}} \sin{\left (\frac{x}{10} \right )} = 0$$

Durch die Multiplikation mit dem Sinus entstehen ja im Prinzip unendlich viele Nullstellen bei den Werten $x_i= +/- \, i \cdot 10 \,\pi \, , \,\, i \in$ ℤ . Die ersten beiden positive Nullstellen des Beitrages vom Sinus lauten z.B $x_0=0$ und $x_1= 10 \,\pi \approx 31.4159$. Die weiteren Nullstellen ergeben sich durch die Nullstellen des Termes $\left( x^4 -q\,x^2 - 5 \right)$, den man einfach mittels einer Variablensubstitution bestimmen kann. In Abhängigkeit vom Parameter $q$ erhält man mittels des Befehles 'solve(...)' die folgenden Lösungen:

In [3]:
Loes_Gl_Nullstelle=solve(Gl_Nullstelle,x)
Loes_Gl_Nullstelle
Out[3]:
$$\left [ 0, \quad 10 \pi, \quad - \sqrt{\frac{q}{2} - \frac{\sqrt{q^{2} + 20}}{2}}, \quad \sqrt{\frac{q}{2} - \frac{\sqrt{q^{2} + 20}}{2}}, \quad - \sqrt{\frac{q}{2} + \frac{\sqrt{q^{2} + 20}}{2}}, \quad \sqrt{\frac{q}{2} + \frac{\sqrt{q^{2} + 20}}{2}}\right ]$$

Für $q=9$ erhalten wir vier reellwertige Lösungen $x_0=0$, $x_1=31.4159265358979$ $x_2=-3.08624979717463$ und $x_3=3.08624979717463$

In [4]:
set_q = 9
for L in Loes_Gl_Nullstelle:
    print(L.subs(q,set_q).evalf())

0
31.4159265358979
-0.724525921248126*I
0.724525921248126*I
-3.08624979717463
3.08624979717463

Drei dieser Nullstellen erkennt man auch in der folgenden Grafik:

In [13]:
plot(f.subs(q,set_q),xlim=(-7,7),ylim=(-2.5,2.5));

Darstellung von $f(x)$, $f^\prime(x)$ und $f^{\prime\prime}(x)$ mittels Jupyter

Stellen Sie nun die Funktion $f(x)$ und die erste und zweite Ableitung der Funktion im Teilintervall $x \in [-7,7]$ grafisch in drei x-y-Diagrammen dar (benutzen Sie dafür wieder den festen Parameterwert $q=9$). Versuchen Sie die Maxima der Funktion (mit $q=9$) mittels Sympy zu bestimmen.

Die erste Ableitung lautet:

In [6]:
diff(f,x)
Out[6]:
$$- \frac{x \left(- q x^{2} + x^{4} - 5\right) e^{- \frac{x^{2}}{4}} \sin{\left (\frac{x}{10} \right )}}{2} + \left(- 2 q x + 4 x^{3}\right) e^{- \frac{x^{2}}{4}} \sin{\left (\frac{x}{10} \right )} + \frac{\left(- q x^{2} + x^{4} - 5\right) e^{- \frac{x^{2}}{4}} \cos{\left (\frac{x}{10} \right )}}{10}$$

bzw. mittels 'simplify()'

In [7]:
diff(f,x).simplify()
Out[7]:
$$\frac{\left(- 20 x \left(q - 2 x^{2}\right) \sin{\left (\frac{x}{10} \right )} + 5 x \left(q x^{2} - x^{4} + 5\right) \sin{\left (\frac{x}{10} \right )} + \left(- q x^{2} + x^{4} - 5\right) \cos{\left (\frac{x}{10} \right )}\right) e^{- \frac{x^{2}}{4}}}{10}$$

Wir stellen uns die Ableitung $f^\prime (x)$ für $q=9$ im Teilintervall $x \in [-7,7]$ grafisch dar.

In [8]:
plot(diff(f,x).subs(q,set_q),xlim=(-7,7),ylim=(-2.5,2.5),ylabel=r"$f^\prime (x)$");

Die zweite Ableitung $f^{\prime \prime} (x)$ lautet:

In [9]:
diff(f,x,x)
Out[9]:
$$\left(2 x^{2} \left(q - 2 x^{2}\right) \sin{\left (\frac{x}{10} \right )} - \frac{x^{2} \left(q x^{2} - x^{4} + 5\right) \sin{\left (\frac{x}{10} \right )}}{4} - \frac{2 x \left(q - 2 x^{2}\right) \cos{\left (\frac{x}{10} \right )}}{5} + \frac{x \left(q x^{2} - x^{4} + 5\right) \cos{\left (\frac{x}{10} \right )}}{10} - 2 \left(q - 6 x^{2}\right) \sin{\left (\frac{x}{10} \right )} + \frac{51 \left(q x^{2} - x^{4} + 5\right) \sin{\left (\frac{x}{10} \right )}}{100}\right) e^{- \frac{x^{2}}{4}}$$

Wir stellen uns die zweite Ableitung $f^{\prime \prime} (x)$ für $q=9$ im Teilintervall $x \in [-7,7]$ grafisch dar.

In [10]:
plot(diff(f,x,x).subs(q,set_q),xlim=(-7,7),ylim=(-3.5,3.5),ylabel=r"$f^{\prime \prime} (x)$");

Wir versuchen nun die Maxima der Funktion (mit $q=9$) mittels Sympy zu bestimmen. Dazu definieren wir zunächst die Gleichung $f^{\prime \prime} (x) = 0$ und benutzen dann, wie sonst auch, den Befehl 'solve(...)':

In [11]:
Gl_Extrema = Eq(diff(f,x).subs(q,set_q),0)
Gl_Extrema
Out[11]:
$$- \frac{x \left(x^{4} - 9 x^{2} - 5\right) e^{- \frac{x^{2}}{4}} \sin{\left (\frac{x}{10} \right )}}{2} + \left(4 x^{3} - 18 x\right) e^{- \frac{x^{2}}{4}} \sin{\left (\frac{x}{10} \right )} + \frac{\left(x^{4} - 9 x^{2} - 5\right) e^{- \frac{x^{2}}{4}} \cos{\left (\frac{x}{10} \right )}}{10} = 0$$
In [12]:
solve(Gl_Extrema,x)

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-12-2ecd30c945e6> in <module>
----> 1 solve(Gl_Extrema,x)

~/anaconda3/lib/python3.7/site-packages/sympy/solvers/solvers.py in solve(f, *symbols, **flags)
   1160     ###########################################################################
   1161     if bare_f:
-> 1162         solution = _solve(f[0], *symbols, **flags)
   1163     else:
   1164         solution = _solve_system(f, symbols, **flags)

~/anaconda3/lib/python3.7/site-packages/sympy/solvers/solvers.py in _solve(f, *symbols, **flags)
   1569                         # solutions in the un-rewritten form below
   1570                         flags['check'] = False
-> 1571                         result = _solve(newf, symbol, **flags)
   1572                         flags['check'] = check
   1573 

~/anaconda3/lib/python3.7/site-packages/sympy/solvers/solvers.py in _solve(f, *symbols, **flags)
   1733 
   1734     if result is False:
-> 1735         raise NotImplementedError('\n'.join([msg, not_impl_msg % f]))
   1736 
   1737     if flags.get('simplify', True):

NotImplementedError: multiple generators [x, tan(x/20)]
No algorithms are implemented to solve equation -10*x*(x**4 - 9*x**2 - 5)*tan(x/20)/(tan(x/20)**2 + 1) + 20*(4*x**3 - 18*x)*tan(x/20)/(tan(x/20)**2 + 1) + (-tan(x/20)**2 + 1)*(x**4 - 9*x**2 - 5)/(tan(x/20)**2 + 1)

... leider ist die Gleichung zu kompliziert um sie auf symbolischem, analytischem Weg zu lösen ...

C++ Programm zur Bestimmung der Hochpunkte der Funktion

Die Berechnung der x-Werte der Hochpunkte der Funktion $f(x)$ wird unter folgendem Link mittels der Methode der Bisektion (Intervallhalbierungsverfahren) in einem C++ Programm berechnet (siehe Musterlösung zum Übungsblatt 7 und das C++ Programm A7_1_Bisektion.cpp).

Es berechneten sich die folgenden beiden x-Werte für die Hochpunkte: $x_{H_1} = -1.8178520203$ und $x_{H_2} = 3.9772608131$. Die berechneten Werte stimmen gut mit den Hochpunkten der Funktion überein, wie man in der folgenden Abbildung sieht.

In [14]:
plot(f.subs(q,set_q),xlim=(-2,4.5),ylim=(-2.5,2.5));
In [ ]: