Demo 6: Ekstremumundersøgelse og optimering

Demo af Christian Mikkelstrup, Hans Henrik Hermansen, Jakob Lemvig, Karl Johan Måstrup Kristensen og Magnus Troen. Revideret marts 2026 af shsp.

from sympy import *
from dtumathtools import *
init_printing()

Til undersøgelsen af extrema for en funktion \(f:\Omega \to \mathbb{R}\) fortæller lærebogen, at de kan findes tre steder:

1. I stationære punkter, altså punkter \(\boldsymbol{x}_0\) hvor \(\nabla f(\boldsymbol{x}_0) = \mathbf 0\)

2. På randen af definitionsmængden, så et punkt \(\boldsymbol{x}_0 \in \partial \Omega \cap \Omega\)

3. I undtagelsespunkter, hvilket vil sige punkter \(\boldsymbol{x}_0\) hvori \(f\) ikke er differentiabel

Karakterisering af lokale ekstrema

Når de stationære punkter for en funktion er fundet, kan de jf. lærebogen karakteriseres ud fra Hessematricens egenværdier i punktet:

• Hvis alle egenværdier er positive, er det stationære punkt et egentligt (lokalt) minimum.

• Hvis de alle er negative, er det et egentligt (lokalt) maksimum.

• Hvis nogle egenværdier har forskelligt fortegn, er det stationære punkt et saddelpunkt, og dermed ikke et extremum.

Bemærk, at hvis en egenværdi er nul i et stationært punkt (mens alle andre egenværdier har samme fortegn), så fortæller denne metode os intet om punktet, og vi vil skulle undersøge punktet på anden vis.

Når karakterisering er nemt (en pæn Hessematrix)

Lad os finde ekstrema for en funktion \(f: \mathbb{R}^2 \to \mathbb{R}\) givet ved:

\[\begin{equation*} f(x,y) = x^3 - 3x^2 + y^3 - 3y^2. \end{equation*}\]

x,y = symbols("x y", real=True)
def f(x,y):
    return x**3 - 3 * x**2 + y**3 - 3 * y**2

f(x,y)

Vi bemærker, at \(f\) er defineret og differentiabel overalt på \(\mathbb{R}^2\). Der er derfor ingen randpunkter eller undtagelsespunkter at undersøge, og ethvert (lokalt) ekstrema vil være at finde i et stationært punkt.

Lad os finde alle stationære punkter:

lign1 = Eq(f(x,y).diff(x),0)
lign2 = Eq(f(x,y).diff(y),0)
lign1,lign2

losn = nonlinsolve([lign1, lign2], [x,y])
losn

og herefter Hessematricen:

H = dtutools.hessian(f(x,y))
H

\[\begin{split}\displaystyle \left[\begin{matrix}6 x - 6 & 0\\0 & 6 y - 6\end{matrix}\right]\end{split}\]

som evalueres i hvert stationært punkt:

[{(x0, y0): H.subs([(x,x0),(y,y0)])} for (x0,y0) in losn]

\[\begin{split}\displaystyle \left[ \left\{ \left( 0, \ 0\right) : \left[\begin{matrix}-6 & 0\\0 & -6\end{matrix}\right]\right\}, \ \left\{ \left( 0, \ 2\right) : \left[\begin{matrix}-6 & 0\\0 & 6\end{matrix}\right]\right\}, \ \left\{ \left( 2, \ 0\right) : \left[\begin{matrix}6 & 0\\0 & -6\end{matrix}\right]\right\}, \ \left\{ \left( 2, \ 2\right) : \left[\begin{matrix}6 & 0\\0 & 6\end{matrix}\right]\right\}\right]\end{split}\]

Da alle disse er diagonalmatricer, læses egenværdierne direkte af deres diagonaler. Vi konkluderer:

• \((0,0)\) er et egentligt (lokalt) maksimum (da vi ser to negative egenværdier).

• \((2,2)\) er et egentligt (lokalt) minimum (da vi ser to positive egenværdier).

• \((0,2)\) og \((2,0)\) er ikke ekstrema men saddelpunkter (da vi ser egenværdier med forskellige fortegn).

Et plot med alle fire stationære punkter:

# Tegn punkterne
liste_med_stat_pt = [Matrix([x0,y0,f(x0,y0)]) for (x0,y0) in losn]
punkter = dtuplot.scatter(liste_med_stat_pt, show=False, rendering_kw={"s" : 100,"color":"red"})

# Tegn fladen f
pf = dtuplot.plot3d(f(x,y),(x,-0.8,2.8),(y,-0.8,2.8),use_cm=True, colorbar=False,show=False,wireframe=True,rendering_kw={"alpha":0.6})
pf.camera = {"azim" : -50, "elev" : 30}

# Samlet plot
pf.extend(punkter)
pf.show()

Ingen globale ekstrema

Bemærk, at vi kun undersøgte lokale extrema i ovenstående. Kunne vi være så heldige, at én af dem også et globalt maksimum eller globalt minimum?

Det ville kræve, at ingen andre punkter giver en større, hhv. lavere, funktionsværdi. Men fordi \(f(x,y)\to\infty\) for fx \(x\to \infty\), da \(x^3\)-leddet gror hurtigere end det negative \(x^2\)-led, og da \(f(x,y)\to-\infty\) for fx \(x\to -\infty\), så er billedmængden af funktionen \(\operatorname{im}(f)=]-\infty,\infty[\,\), og funktionen har dermed ingen globale ekstrema.

Når karakterisering ikke er så nemt

Lad os nu undersøge ekstrema for en anden funktion \(f:\mathbb{R}^2 \to \mathbb{R}\) givet ved

\[\begin{equation*} f(x,y) = x^4 + 4x^2y^2 + y^4 -4x^3 - 4y^3 + 2. \end{equation*}\]

x,y = symbols('x y', real=True)
def f(x, y):
    return x**4 + 4*x**2*y**2 + y**4 - 4*x**3 - 4*y**3 + 2

f(x,y)

Igen er definitionsmængden hele \(\mathbb{R}^2\), hvorpå \(f\) er differentiabel overalt, så igen er der ingen rand ej heller undtagelsespunkter at undersøge, og vi fokuserer kun på stationære punkter.

Alle stationære punkter er:

lign1 = Eq(f(x,y).diff(x),0)
lign2 = Eq(f(x,y).diff(y),0)
lign1,lign2

losn = nonlinsolve([lign1,lign2], [x, y])
losn

Lad os få disse udskrevet som decimaltal:

[(N(xlosn),N(ylosn)) for (xlosn,ylosn) in losn]

Vi ser ikke-reelle løsninger, selvom vi tidligere afgrænsede vores variabler til de reelle tal med real=True. Den ærgerlige årsag er, at nonlinsolve() ignorerer, hvordan variablerne er defineret. Vi er derfor nødt til manuelt at udtrække de løsninger fra løsningsmængden, der ligger inden for definitionsmængden. Der er kun fire reelle løsninger iblandt dem:

stat_pt = list(losn)[:-2] # Udelad de sidste to elementer (de ikke-reelle løsninger) med "list slicing"
stat_pt

liste_med_stat_pt = [Matrix([x0,y0,f(x0,y0)]) for (x0,y0) in stat_pt]
punkter = dtuplot.scatter(liste_med_stat_pt, show=False, rendering_kw={"color":"red", "s" : 30})

pf = dtuplot.plot3d(f(x,y),(x,-2, 5),(y,-2, 4),xlim=(-2,5),ylim=(-2,5),zlim=(-30,15), show=False, rendering_kw={"alpha":0.5})
pf.camera = {"azim" : -161, "elev" : 5}
(pf + punkter).show()

Hessematricen

H = dtutools.hessian(f(x,y))
H

\[\begin{split}\displaystyle \left[\begin{matrix}12 x^{2} - 24 x + 8 y^{2} & 16 x y\\16 x y & 8 x^{2} + 12 y^{2} - 24 y\end{matrix}\right]\end{split}\]

De stationære punkter indsættes én ad gangen:

Hesse_matricer = [H.subs([(x,x0),(y,y0)]) for x0,y0 in stat_pt]
Hesse_matricer

\[\begin{split}\displaystyle \left[ \left[\begin{matrix}0 & 0\\0 & 0\end{matrix}\right], \ \left[\begin{matrix}72 & 0\\0 & 36\end{matrix}\right], \ \left[\begin{matrix}-4 & 16\\16 & -4\end{matrix}\right], \ \left[\begin{matrix}36 & 0\\0 & 72\end{matrix}\right]\right]\end{split}\]

Ikke alle Hessematricer er diagonale, så vi beregner deres egenværdier med eigenvals:

Egenv_Hesse_matricer = [h.eigenvals() for h in Hesse_matricer]

display(*zip(stat_pt,Egenv_Hesse_matricer))

Heraf aflæser vi:

• \((0,3)\) og \((3,0)\) er begge egentlige (lokale) minima (positive egenværdier, \(\lambda=36,72\))

• \((1,1)\) er et saddelpunkt (egenværdier med modsatte fortegn, \(\lambda=-20,12\))

Det sidste stationære punkt \((0,0)\) kan desværre ikke karakteriseres med denne metode, da vi her har en egenværdi på \(0\). Dette punkt skal undersøges særskildt på anden vis.

Særundersøgelse af problematisk stationært punkt

En metode til særundersøgelse af punktet, er at nærme sig det fra forskellige retninger på funktionen. Dette kan gøres ved, at der laves en restriktion af \(f\) langs forskellige kurver i definitionsmængden. Hvis \((0,0)\) vitterligt er et ekstremumpunkt for \(f\), så vil det være et ekstremumpunkt for enhver kurve på \(f\), der går igennem punktet (er det fx et maksimum for \(f\), så er det et maksimum for alle kurver igennem det). Skulle vi være heldige at finde blot ét modeksempel, så kan vi straks konkludere, at \(f\) ikke har ekstremum i punktet.

Lad os fx undersøge, hvordan \(f\) opfører sig langs linjen \(y=x\) (altså, vi undersøger restriktionen \(f(x,x)\)):

dtuplot.plot(f(x,x),f(0,0),(x,-0.5,1.2),axis_center="auto")

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f91d8bdf9d0>

Vi har tilføjet en vandret linje for at angive funktionsværdien i punktet, \(f(0,0)\). Husk på, at hver \(x\)-værdi på dette plot repræsenterer et \((x,x)\)-punkt på \(f\); det er altså i \(x=0\) på denne graf, at vi finder \(f\)s stationære punkt \((0,0)\). Visuelt ligner dette punkt et saddelpunkt langs \(y=x\). Lad os undersøge forskellen \(f(x,x) - f(0,0)\), når vi nærmer os punktet fra hver side langs \(y=x\):

forskel = f(x,x) - f(0,0)
forskel.factor()

Vi ser, at

\[f(x,x) - f(0,0) > 0 \text{ for } x < 0\,, \quad\text{og}\quad f(x,x) - f(0,0) < 0 \text{ for } 0 < x < \frac{4}{3}.\]

Da \(f(x,x)\) antager værdier, der både er større og mindre end \(f(0,0)\) for \((x,x)\) vilkårligt tæt på \((0,0)\), så kan \(f\) ikke have ekstremum i \((0,0)\). \((0,0)\) er et saddelpunkt på linjen \(y=x\) og dermed også for \(f\).

Globale ekstrema

Lad os nu betragte en funktion \(f:\Omega \to \mathbb{R}\) med forskriften

\[\begin{equation*} f(x,y) = 3 + x - x ^ 2 - y ^ 2. \end{equation*}\]

def f(x,y):
    return 3 + x - x ** 2 - y ** 2

f(x,y)

\(f\) er veldefineret på hele \(\mathbb{R}^2\), men denne gang begrænser vi definitionsmængden ved

\[\begin{equation*} \Omega = \Bigl\{(x,y) \in [0,2]\times[-1,1]\quad |\quad y \leq 1-x\Bigr\}, \end{equation*}\]

altså som mængden, der omkranses af linjerne \(y = -1\), \(y = 1-x\) og \(x=0\).

M = dtuplot.plot_implicit(Eq(x,0),Eq(y,1-x),Eq(y,-1), (1-x >= y) & (y >= -1) & (x >= 0),(x,-0.1,2.1),(y,-1.1,1.1), show=False, adaptive=True)
M.legend = False
M.show()

\(f\) er differentiabel på hele \(\Omega\), så der er ingen undtagelsespunkter. Derfor er den også kontinuert på hele \(\Omega\). Da \(\Omega\) er afsluttet og begrænset, har \(f\) et globalt minimum og et globalt maksimum på \(\Omega\). Disse to værdier antages enten i stationære punkter eller på randen.

Da \(\Omega\) er en (kurve)sammenhængende mængde, ved vi, at billedmængden/værdimængden er \(f(\Omega) = [m;M]\), hvor \(m\) er det globale minimum og \(M\) det globale maksimum.

Her et hurtigt blik på \(f\)s niveaukurver:

niveaukurve = dtuplot.plot_contour(f(x,y),(x,-0.1,2.1),(y,-1.1,1.1),show=False)
niveaukurve.extend(dtuplot.plot_implicit(Eq(x,0),Eq(y,1-x),Eq(y,-1),(x,-0.1,2.1),(y,-1.1,1.1),show=False))
niveaukurve.show()

Undersøgelse af stationære punkter

Vi sætter de partielt afledte lig med nul og løser systemet:

lign1 = Eq(f(x,y).diff(x),0)
lign2 = Eq(f(x,y).diff(y),0)
losn = nonlinsolve([lign1,lign2], [x,y])
losn

Her er kun én løsning, og det er et punkt inden for definitionsmængden \(\Omega\). Derfor har \(f\) ét stationært punkt, hvilket er en kandidat til at være globalt ekstremum. Vores søgning er ikke slut endnu, så lad os placere dette punkt i en kandidatliste, der kan opdateres løbende, som yderligere kandidater dukker op:

kandidater = set(losn)
kandidater

Randundersøgelse

Næste trin er at finde mulige kandidater for globalt ekstremum på randen af definitionsmængden. Grundet hjørnerne, er randen kun stykkevist differentiabel. Vi undersøger hvert differentiabelt randstykke individuelt.

For at gøre det, laver vi en restriktion af \(f\) til hvert randstykke, ét ad gangen, således at vi kan undersøge stykkerne for stationære punkter. De tre restriktioner dannes ved \(f(0,y)\), \(f(x,1-x)\) og \(f(x,-1)\). Lad os lægge ud med at plotte dem:

dtuplot.plot(f(x,-1),(x,0,2), title="Randrestriktion f(x, -1)")

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f91d8a49610>

dtuplot.plot(f(x,1-x),(x,0,2), title="Randrestriktion f(x, 1-x)")

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f91d88f7910>

dtuplot.plot(f(0,y), (y,-1,1), ylim = (0,3),aspect="equal", title="Randrestriktion f(0, y)")

<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f91d8563f10>

Plottet tyder på et stationært punkt på hvert randstykke. For at være sikre differentierer vi dem og sætter resultaterne lig nul:

lodretStykke = solve(f(0,y).diff(y))
lodretStykke

vandretStykke = solve(f(x,-1).diff(x))
vandretStykke

vinkletStykke = solve(f(x,1-x).diff(x))
vinkletStykke

Husk på, at når vi løser et udtryk uden eksplicit at angive højresiden eller variablen, så antager Sympy, at højresiden skal være nul, og at variablen skal være den eneste tilgængelige variabel. Med solve(f(0,y).diff(y)) løses altså ligningen \(\frac{\mathrm df(0,y)}{\mathrm dy}=0\) for \(y\), og vi ser herover, at resultatet blev \(y=0\).

Vi har ganske rigtigt fundet én løsning for hvert randstykke, som alle befinder sig på deres respektive randstykke (altså, inden for definitionsmængden af randstykket) og derfor inden for \(f\)s definitionsmængde. De er derfor alle tre stationære punkter på randen af \(f\) og mulige kandidater til globalt ekstremum.

Lad os rekonstruere dem som punkter i \(\Omega\subset \Bbb R^2\) og tilføje dem til vores kandidatliste:

kandidater.update(set([(0,y0) for y0 in lodretStykke]))
kandidater.update(set([(x0,-1) for x0 in vandretStykke]))
kandidater.update(set([(x0,1-x0) for x0 in vinkletStykke]))
kandidater

Glem ikke hjørnerne

Har vi nu undersøgt hele definitionsmængden \(\Omega\) for \(f\)? Har vi betragtet alle punkter?

Nej! Vi mangler hjørnerne på \(\Omega\), som er \((0,1)\), \((0,-1)\) og \((2,-1)\). Disse punkter ligger på randen, men da randrestriktionen af \(f\) ikke er differentiabel dér, blev de ekskluderet i vores randundersøgelse ovenfor, da hvert randstykke blev betragtet som en åben mængde og dermed uden endepunkter.

Heldigvis er enkelte, manglende punkter nemme at håndtere. Tilføj dem bare til kandidatlisten. Det er blot vigtigt, at vi ikke glemmer dem.

kandidater = list(kandidater) + [(0,1),(0,-1),(2,-1)]
kandidater

Den afsluttende sammenligning

Vi har nu samlet en lang liste af mulige kandidater til globale ekstrema efter en grundig, udtømmende søgning gennem definitionsmængden. Det sidste, der mangler, er en sammenligning af deres funktionsværdier - enten ved manuelt at taste punkterne fra kandidatlisten ind:

f(3/4,1/4),f(1/2,-1),f(1/2,0),f(0,0),f(0,1),f(0,-1),f(2,-1)

eller, hvis du vil undgå risikoen for tastefejl, at lade Sympy udtrække og evaluere alle punkterne i én kodelinje:

f_vaerdier = [f(x0,y0) for x0,y0 in kandidater]
f_vaerdier

Da vi har udtømt definitionsmængden, kan vi være sikre på, at det punkt i vores kandidatliste, der giver den største funktionsværdi, må være funktionens globale maksimum, mens kandidaten, der giver den mindste funktionsværdi, må være dens globale minimum. Vi kan direkte se i listen, hvilke der er tale om, eller vi kan lade Sympy udvælge den største og mindste funktionsværdi:

minimum = min(f_vaerdier)
maksimum = max(f_vaerdier)
minimum, kandidater[f_vaerdier.index(minimum)], maksimum, kandidater[f_vaerdier.index(maksimum)]

Altså:

• Det globale minimum for \(f\) er \(0\) og antages i punktet \((2,-1)\) (et af hjørnerne)

• Det globale maksimum for \(f\) er \(\frac{13}{4}\) og antages i punktet \((\frac{1}{2},0)\) (det indre stationære punkt)

Som samlet konklusion har vi, at billedmængden af \(f\) er \(\mathrm{im} f=[0,\frac{13}{4}]\).

Et plot af grafen for \(f\) med definitionsmængden markeret samt de to globale ekstrema indtegnet:

u = symbols("u")
pf = dtuplot.plot3d(f(x,y),(x,0,2),(y,-1,1),show=False,rendering_kw={"alpha":0.7})
punkter = dtuplot.scatter([Matrix([2,-1,0]),Matrix([1/2,0,13/4])],show=False,rendering_kw={"color":"red","s":20})
l1 = dtuplot.plot3d_parametric_line(u,-1,f(u,-1),(u,0,2),use_cm=False,show=False,rendering_kw={"color":"red","linewidth":2})
l2 = dtuplot.plot3d_parametric_line(0,u,f(0,u),(u,-1,1),use_cm=False,show=False,rendering_kw={"color":"red","linewidth":2})
l3 = dtuplot.plot3d_parametric_line(u,1-u,f(u,1-u),(u,0,2),use_cm=False,show=False,rendering_kw={"color":"red","linewidth":2})
kombineret = (pf + punkter + l1 + l2 + l3)
kombineret.camera = {"azim" : 148, "elev" : 40}
kombineret.legend = False

kombineret.show()