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:
I stationære punkter, altså punkter \(\boldsymbol{x}_0\) hvor \(\nabla f(\boldsymbol{x}_0) = \mathbf 0\)
På randen af definitionsmængden, så et punkt \(\boldsymbol{x}_0 \in \partial \Omega \cap \Omega\)
I undtagelsespunkter, hvilket vil sige punkter \(\boldsymbol{x}_0\) hvori \(f\) ikke er differentiabel
Karakterisering af lokale ekstrema#
Stationære punkter er løsninger til ligningen \(\nabla f(\boldsymbol x)=\mathbf 0\). Dette er egentlig et ligningssystem, hvor hver partielt afledte er sat lig nul. Når det er løst, og funktionens stationære punkter er fundet, har vi ofte behov for at afgøre, hvilken type de hver især er - vi skal karakterisere de stationære punkter. Til dette formål henviser lærebogen os til at betragte Hessematricens egenværdier i hvert stationært punkt:
Hvis alle egenværdier er positive, er det stationære punkt et egentligt (lokal) minimum.
Hvis de alle er negative, er det et egentligt (lokal) 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ære punkt, 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:
x,y = symbols("x y", real=True)
f = x**3 - 3 * x**2 + y**3 - 3 * y**2
Vi bemærker, at \(f\) er defineret og differentiabel overalt i hele \(\mathbb{R}^2\). Der er derfor ingen randpunkter eller undtagelsespunkter at undersøge. Ethvert (lokalt) ekstrema vil derfor være at finde i et stationært punkt.
Lad os finde alle stationære punkter:
lign1 = Eq(f.diff(x),0)
lign2 = Eq(f.diff(y),0)
display(lign1,lign2)
losn = nonlinsolve([lign1, lign2], [x,y])
losn
Og herefter Hessematricen:
H = dtutools.hessian(f)
H
som evalueres i hvert stationært punkt:
[{(x0, y0): H.subs([(x,x0),(y,y0)])} for (x0,y0) in losn]
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 (da vi ser egenværdier med forskellige fortegn). Disse er derimod saddelpunkter.
Et plot med alle fire stationære punkter:
# Tegn punkterne
liste_med_stat_pt = [Matrix([x0,y0,f.subs([(x,x0),(y,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,-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
x,y = symbols('x y', real=True)
f = x ** 4 + 4 * x**2 * y ** 2 + y ** 4 - 4 * x ** 3 - 4 * y ** 3 + 2
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.
Vi finder alle stationære punkter:
lign1 = Eq(f.diff(x),0)
lign2 = Eq(f.diff(y),0)
lign1,lign2
lign = nonlinsolve([lign1,lign2], [x, y])
lign
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 uheldige årsag er, at nonlinsolve() ignorerer, hvordan variablerne er defineret. Vi skal derfor manuelt udtrække de løsninger, der ligger inden for vores ønskede definitionsmængde, fra løsningsmæ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.subs([(x,x0),(y,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,-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)
H
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
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#
Først og fremmest er funktionsværdien i dette problematiske stationære punkt \((0,0)\),
f.subs({x: 0, y: 0})
lig med \(f(0,0)=2.\)
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.subs(y,x),2,(x,-0.5,1.2),axis_center="auto")
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b79b58dd0>
Plottet har her fået tilføjet en vandret linje, der repræsenterer funktionsværdien i punktet. 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)\), som vi nærmer os punktet fra hver side langs \(y=x\):
tmp = f.subs(y,x) - f.subs([(x,0),(y,0)])
tmp.factor()
Vi ser at
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
\(f\) er veldefineret på hele \(\mathbb{R}^2\), men denne gang begrænser vi definitionsmængden ved
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:
f = 3 + x - x ** 2 - y ** 2
niveaukurve = dtuplot.plot_contour(f,(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.diff(x),0)
lign2 = Eq(f.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.subs(y, -1),(x,0,2), title="Randrestriktion f(x, -1)")
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b7772d010>
dtuplot.plot(f.subs(y, 1-x),(x,0,2), title="Randrestriktion f(x, 1-x)")
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b5697fed0>
dtuplot.plot(f.subs(x,0), (y,-1,1), ylim = (0,3),aspect="equal", title="Randrestriktion f(0, y)")
<spb.backends.matplotlib.matplotlib.MatplotlibBackend at 0x7f0b56ab9350>
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.subs(x,0).diff(y))
lodretStykke
vandretStykke = solve(f.subs(y,-1).diff(x))
vandretStykke
vinkletStykke = solve(f.subs(y,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.subs(x,0).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 er 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:
f_vaerdier = [f.subs([(x,x0),(y,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 trække punkterne ud fra listen manuelt for derefter at sætte dem ind i \(f\) eller lade Sympy gøre det for os:
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 \([0,\frac{13}{4}]\).
Et plot med grafen for \(f\) men definitionsmængden markeret samt de to globale ekstrema angivet:
u = symbols("u")
pf = dtuplot.plot3d(f,(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.subs({x:u,y:-1}),(u,0,2),use_cm=False,show=False,rendering_kw={"color":"red","linewidth":2})
l2 = dtuplot.plot3d_parametric_line(0,u,f.subs({x:0,y:u}),(u,-1,1),use_cm=False,show=False,rendering_kw={"color":"red","linewidth":2})
l3 = dtuplot.plot3d_parametric_line(u,1-u,f.subs({x:u,y: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()
