Uge 1: Funktioner, Partielle afledede, Gradient-vektoren#

Demo af Christian Mikkelstrup og Hans Henrik Hermansen

from sympy import *
from dtumathtools import *
init_printing()

Velkommen tilbage efter jul og januar, og velkommen til foråret i mat1. Der kommer til at være en helt masse nyt pensum, og blandt andet en del 3D-plots! Til dette har vi udviklet dtumathtools, som vil følge jer i løbet af foråret. Den indeholder dtuplot som skal bruges til at plotte, samt flere gode hjælpefunktioner. I burde have dtumathtools instaleret fra mat1a, hvis ikke så kør kommandoen

pip install dtumathtools 

i terminalen.

Efter dette kan resten af demoen nydes. Vi starter med noget kendt materiale, nemlig funtioner af en variabel.

Funktioner (af en variabel)#

Man kan definere funktionen \(f: \mathbb{R} \to \mathbb{R}\), \(f(x)=x \mathrm{e}^x\) som en Python-funktion med følgende kendte kommando:

def f(x):  
    return x * exp(x)

Funktionen evalueres i punktet \(x=-2\) med kommandoen:

f(-2)
../_images/5143a366292d0b8d11a3379e75854f4443da9c761b4e42b6a74d180e08bc3d79.png

hvis numeriske værdi er:

f(-2).evalf()
../_images/0366f0d81388fb2f714c7ad32d87eaf184af7bf112e770b5656ef8bb2cd79bc5.png

Det er dog sjældent nødvendigt at definere funktioner ved def-kommandoen, og vi vil ofte blot arbejde direkte med funktions-udtrykket:

x = symbols('x')          # nødvendigt for at kunne bruge x som symbolsk variabel
f_expr = x * exp(x)
f_expr
../_images/bb59df58659ff15d76d3315db722c41be2564d1029f5061b6a15bbc2bc40d045.png

som evalueres med:

f_expr.subs(x, -2)
../_images/5143a366292d0b8d11a3379e75854f4443da9c761b4e42b6a74d180e08bc3d79.png

Funktionen kan fx differentieres ved:

f_maerke = f_expr.diff(x)
f_maerke
../_images/196f67a172cde873c4d8ac79a56d769d3686a88d00f1b8e09d1e1851f0f389f5.png

og vi kan undersøge grænseovergange for \(x \to -\infty\), \(x \to \infty\) og \(x \to -2\) med:

f_expr.limit(x, -oo), f_expr.limit(x, oo), f_expr.limit(x, -2)
../_images/06902cfb75e0c87555f02cb97e701e9b1b30fe56f3415444175816aa8df020a2.png

Da funktionen er kontinuert overalt er det ikke overraskende at \(\lim_{x \to -2} f(x) = f(-2)\).

Funktionen og den afledte plottes ved:

plot(f_expr, f_maerke, (x, -5, 1))
../_images/b938880ec60078719c34ac467eaeab410bdc6a54935396179b0eb458a75df265.png
<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7fa622ba89d0>

En lidt mere kompliceret eksempel er den stykkevis definerede funktion \(g: \mathbb{R} \to \mathbb{R}\):

\[\begin{equation*} g(x) = \begin{cases} -x & x <0 \\ \mathrm{e}^x & x \ge 0 \end{cases} \end{equation*}\]

Som Python-funktion er dette:

def g(x):
    if x < 0:
        return -x
    else:
        return exp(x)

Igen er det dog smartere at bruge funktions-“udtrykket” i Sympy:

g_expr = Piecewise((-x, x < 0), (exp(x), x >= 0))

som igen evalueres ved:

g_expr.subs(x, -2)
../_images/b1fac8f1225f03d1681446e2a59ec935caa3f7b0c81185479a10e26dd6c63c3b.png

og plottes ved:

plot(g_expr,(x,-5,2), ylabel='g(x)')
../_images/b1ddf369032324f09860dcc79665036edc6447851e5dbf42b6720f334f9fff48.png
<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7fa601d08f70>

Vi ser at funktionen ikke er kontinuert i \(x=0\). Bemærk at man ikke kan bruge Python/CAS til at bevise at funktionen er diskontinuert i \(x=0\). Dette kræver et såkaldt epsilon-delta-argument og her er Python/CAS ikke nyttig. Sympy vil fx gerne differentiere funktionen:

g_expr.diff(x)
\[\begin{split}\displaystyle \begin{cases} -1 & \text{for}\: x < 0 \\e^{x} & \text{otherwise} \end{cases}\end{split}\]

men vi bør bemærke at funktionen ikke er differentiabel i \(x=0\) da den ikke er kontinuert i \(x=0\).

Funktioner af flere variable: Partielt afledte ved brug af diff#

For funktioner af flere variable introducerer vi partielle afledte, samt hvordan vi kan benytte dem. For at vise hvordan man kan få de partielt afledte, kan vi kigge på funktionen:

x, y = symbols('x y')
f = x*y**2+x
f
../_images/1826ecb9e7061052d604b7e74e71fce0c2df70b654c6589cce14f54f8fa26712.png

og finde de afledte med kommandoen som vi også brugte sidste semester

f.diff(x), f.diff(y)
../_images/d49ce783966f42544970fd731cfcc6fe6a24935bebcefd9ba2d3856e9ed8f0cd.png

Hver af disse udtryk kan vi jo igen differentiere mht \(x\) og \(y\). Dette giver følgende fire funktioner:

f.diff(x).diff(x), f.diff(x).diff(y), f.diff(y).diff(x), f.diff(y).diff(y)
../_images/c6708f51a80726290ada31bf2c7ee4a8de4f64fdbe1e7f5eda8c7c500b226ad0.png

Dette kaldes de afledte af anden orden og kan udregnes direkte ved:

f.diff(x,2), f.diff(x,y), f.diff(y,x), f.diff(y,2)
../_images/c6708f51a80726290ada31bf2c7ee4a8de4f64fdbe1e7f5eda8c7c500b226ad0.png

Vi kan så indsætte værdier i denne, for eksempel \(\frac{\partial}{\partial x}f(x,y)\) taget i \((-2,3)\) ved

f.diff(x).subs({x:-2,y:3})
../_images/ade62f20a7bcb96c8d3ed92debf0a75ef39435164d2934af8a560d7b030168c3.png

Eller \(\frac{\partial}{\partial x\partial y}f(x,y)\), taget i \((5,-13)\) ved

f.diff(x,y).subs({x:5,y:-13})
../_images/970dc01d48cda10258dc783b17ed4b915b5add84cbc606c08feb206c1895dfb6.png

Plots#

Orientering af grafer#

Vi skal nu, for første gang til at plotte funktioner af flere variable, og altså i 3D! Her er et valg man skal tage, for man har nemlig mulighed for at rotere et plot rundt, så man kan se det fra flere vinkler. Hvis ikke man gør noget, så vælger dtuplot en vinkel for os, men hvis man ønsker at se grafen fra en bestemt vinkel kan camera anvendes. Prøv herunder at ændre værdierne for elev og azim.

f = 4-x**2-y**2

p=dtuplot.plot3d(f, (x,-3,3),(y,-3,3), camera = {"elev": 25, "azim": 45})
../_images/9cac37f48acbd6a9724d0f1c1ca95affb4e9b35c2dc0ed14f37210a13d5bfaf8.png

Ovenstående kommando laver plottet som en statisk PNG-fil, hvilket er smart hvis man skal printe Notebook’en eller eksportere til PDF. Statiske plots fås ved ikke at gøre noget eller ved brug af %matplotlib inline.

Hvis man i stedet kører %matplotlib qt (i den følgende blok udkommenteret, men prøv at fjerne udkommenteringen #), aktiverer man interaktive plots. Alle efterfølgende plots “popper” nu ud af VS Code, hvorefter man man rotere plottet rundt og se det fra flere vinkler! Prøv derefter at plotte 3D plottet igen!

#%matplotlib qt

Om interaktive plots#

Bemærk: %matplotlib qt virker normalt kun hvis man kører fx VS Code på egen laptop. Hvis man fx kører Python på en online server i browseren som Google Colab, vil %matplotlib qt ikke virke. Her kan man prøve widgets i stedet for: %matplotlib ipympl. Det kræver at man lige installerer pakken ipympl. Samlet oversigt:

# Fjern udkommentering for den backend der ønskes brugt
# %matplotlib inline        # statisk plots
# %matplotlib qt            # QT (cute) interaktivt pop-ud plots
# %matplotlib ipympl        # Widget/ipynpl interaktivt inline plots (ikke så stabil som QT og kan kræve restart af kernel)
# %matplotlib --list        # liste over alle backends  

Æstetik#

Når man vil ændre på æstetikken af et plot bruges rendering_kw={....} som argument, og det kan se lidt underligt ud. Dette er blot hvilke æstetiske (rendering) indstillinger der skal bruges, eksempelvis color, alpha, osv. Man kan også i det fleste tilfælde “bare” skrive {....}, og så ved den godt at det er det æstetiske, men det er mere tydeligt at skrive det med.

dtuplot.plot3d(f, (x,-3,3),(y,-3,3), wireframe = True, rendering_kw = {"color": "red", "alpha": 0.5})
../_images/8113b8a4f6737eec695e26c5b7ef61f6c19a2a3fcb4ddf54d5e9ef1a850d2f1c.png
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7fa6018fb820>

Nogle æstetiske valg er dog specielle nok til at få lov til at stå alene, fx wireframe som giver grafen et gitter, eller use_cm som herunder angiver værdien i hvert punkt med en farve.

p=dtuplot.plot3d(f, (x,-3,3),(y,-3,3), use_cm=True, legend=True)
../_images/74560ab1df661a1a9f20c2a63cb64549f88667c76d7c5d87a0256ad6d9ffd6d5.png

Niveaukurver#

Vi kan også plotte højdelinjer, altså et 2D plot over en 3D struktur ved:

dtuplot.plot_contour(f, (x,-3,3),(y,-3,3), is_filled=False)
../_images/767b6f873b6546d2a7542a9c2af9b3f990eb5cd4f8f7ace3b554c6ed8ddf4194.png
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7fa6017a89a0>

Og hvis vi vil bestemme hvilke højder der vises, kan vi bruge,

zvals = [-2,-1,0,1]
dtuplot.plot_contour(f, (x,-3,3),(y,-3,3), rendering_kw={"levels":zvals, "alpha":0.5}, is_filled=False)
../_images/6363470bbb687ba92aaccbe6fe0c883f33a81c0e77a95de8604d9ab6f8a00984.png
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7fa601bf3a60>

Gradientvektorfelter#

Vi kigger nu på funktionen \(f: \mathbb{R}^2 \to \mathbb{R}\):

(5)#\[\begin{equation} f(x,y)=\cos(x)+\sin(y). \end{equation}\]

hvis graf i 3D visualiseres ved:

f = cos(x)+sin(y)
p = dtuplot.plot3d(f, (x,-pi/2,3/2*pi),(y,0,2*pi),use_cm=True, camera={"elev":45, "azim":-65}, legend=True)
../_images/237c3df74a3c6e0a854487607738c11a111828bfc6a9a511dec60b34ba86b709.png

Gradienten for \(f\) taget i punktet \((x,y)\) er en vektor, som har symbolet \(\nabla f(x,y)\) (\(\nabla\) kaldes nabla). Den er sammensat af de to partielt afledte,

nf = Matrix([f.diff(x), f.diff(y)])
nf
\[\begin{split}\displaystyle \left[\begin{matrix}- \sin{\left(x \right)}\\\cos{\left(y \right)}\end{matrix}\right]\end{split}\]

Gradienten kan også fås med dtutools.gradient, bemærk dog at man med denne funktion ikke altid har magt over hvilken rækkefølge de variable tages i.

dtutools.gradient(f)
\[\begin{split}\displaystyle \left[\begin{matrix}- \sin{\left(x \right)}\\\cos{\left(y \right)}\end{matrix}\right]\end{split}\]

Vi kan opfatte gradienten som en funktion \(\nabla f: \mathbb{R}^2 \to \mathbb{R}^2\) (som kaldes et gradientvektorfelt), som let plottes ved:

dtuplot.plot_vector(nf, (x,-pi/2,3/2*pi),(y,0,2*pi),scalar=False)
../_images/4d211a828c7cdc61665191e1c6b759201ccebb85487e4856210687a1be8351ac.png
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7fa5f9d79520>

Eller hvis det skal være lidt flottere (her er rendering_kw splittet, så man kan specificere for pile og contours seperat),

dtuplot.plot_vector(nf, (x,-pi/2,3/2*pi),(y,0,2*pi),
    quiver_kw={"color":"black"},
    contour_kw={"cmap": "Blues_r", "levels": 20},
    grid=False, xlabel="x", ylabel="y",n=15)
../_images/8945a7c59f2c92eabbeb3bf9c054989725e4e572d2c4ea7e7f621e4e954ee5b4.png
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7fa5f94d7e80>