1.6. Olinj¨ ara problem

forts. p˚a f¨oreg˚aende f¨orel¨asning:

Minsta kvadratmetoden kan ocks˚a till¨ampas p˚a funktioner som inte har ett linj¨art beroende av parametrarna.

Vi skall b¨orja med att beskriva i korthet den statistiska princip, som ligger till den grund f¨or den tidigare beskrivna minsta kvadratmetoden.

Antag, att vi ¨onskar best¨amma parametervektorn x = {xi, i = 1, 2, . . . , m} i modellekvationen y(t) = f (t, x) utg˚aende fr˚an observationerna

y_i = y(t_i) + e_i (i = 1, 2, . . . , n) (n > m),

d¨ar m¨atfelen ei antas vara oberoende, normalf¨ordelade slumptal med medelv¨ardet 0 och standardavvikelsen σ_i.

Detta inneb¨ar, att sannolikheten f¨or m¨atv¨ardet yi har normalf¨ordelningen

P (y_i) = 1 σ_i√

2πexp (

−¹₂[y_i − f (t_i, x)]² σ_i²

) .

Om m¨atv¨ardena ¨ar oberoende, s˚a kommer sannolikheten f¨or att erh˚alla en best¨amd serie observationsdata y_i, i = 1, 2, . . . , n att vara produkten av de enskilda sannolikheterna

P_y = P (y₁)P (y₂) . . . .

P_y ¨ar givetvis i allm¨anhet en funktion av parametrarna xi (i = 1, 2, . . . , m). Genom att inf¨ora beteck- ningen

χ²(x) =

Xn i=1

1

σ_i²[y_i − f (ti, x)]², s˚a finner man, att den totala sannolikheten kan uttryckas

P_y(x) = 1

σ₁σ₂ . . . σ_n(2π)^n/2e^−χ2/2.

Enligt maximeringsprincipen (se f¨orel¨asn. 10 i kursens f¨orsta del) b¨or parametrarna v¨aljas s˚a, att P_y blir s˚a stor som m¨ojligt. Detta leder till villkoret χ² = minimum, vilket ¨ar ekvivalent med minsta kvadratmetoden¹. Observera, att denna h¨arledning bygger p˚a, att m¨atv¨ardena ¨ar normalf¨ordelade, vilket inte alltid beh¨over vara fallet.

1C.F. Gauss: Theoria Motus..., Art. 179

Som ett uttryck f¨or anpassningens godhet anv¨ander man ofta testfunktionen

S = Xn

k=1

w_ke²_k = Xn

k=1

w_k[y_k − f (tk, x)]²,

ist¨allet f¨or χ². w_k betecknar m¨atningarnas vikter (wk ∝ ¹

σ2k

).

I matrisform kan denna funktion uttryckas S = e⁰W e, d¨ar W ¨ar en diagonal viktsmatris (Wij = w_iδ_ij), och e en felvektor (e_i = y_i − f (t_i, x)).

Om ifr˚agavarande minsta kvadratproblem ¨ar linj¨art, s˚a kan denna felvektor uttryckas i formen e = y −Ax, d¨ar A ¨ar en n × m–matris med konstanta koefficienter (observera, att n > m).

I minimet g¨aller ^∂S

∂xk = 0. Genom att till¨ampa detta villkor p˚a testfunktionen S f˚as

0 = ∂S

∂x_k = −2 Xn

i=1

w_i(y_i − Xm

j=1

A_ijx_j)A_ik,

som leder till ekvationssystemet Xn

i=1

Xm j=1

A_ikw_iA_ijx_j = Xn

i=1

A_ikw_iy_i (k = 1, 2, . . . , n)

(Gauss’ normalekvationer).

Dessa ekvationer kan ocks˚a framst¨allas i matrisform

A⁰W Ax = A⁰W y.

Om matrisen A⁰W A (som ¨ar en symmetrisk m × m–matris) inte ¨ar singul¨ar, s˚a ¨ar x₀ = (A⁰W A)⁻¹A⁰W y

en station¨ar punkt av S. Denna matris ¨ar relaterad till Hesses matris (efter Otto Hesse, 1811-1874) f¨or S, som ¨ar H = 2A⁰W A. Man kan visa, att x₀ ¨ar ett minimist¨alle, ifall H ¨ar positivt definit. Detta ¨ar analogt med det endimensionella fallet d¨ar x0 ¨ar ett minimist¨alle, om f⁰(x₀) = 0 och f⁰⁰(x₀) > 0.

Ifall f beror olinj¨art av x, s˚a kan minsta kvadratproblemet (i princip) l¨osas genom successiva approxima- tioner, d¨ar man vid varje steg l¨oser ett linj¨art minsta kvadratproblem.

1.6. Olinj¨ ara problem

Det f¨orekommer ofta, att man vill best¨amma v¨ardet av parametrar, som p˚a ett olinj¨art s¨att beror av m¨atbara storheter. Vi har t.ex. en kemisk reaktion, d¨ar jonkoncentrationen vid en viss tidpunkt t antar ett v¨arde som kan uttryckas i formen f (t) = 10e^−3t + 2e^−5t. Antag, att vi vill veta n¨ar jonkoncentrationen nedg˚att till h¨alften av sitt ursprungliga v¨arde. Emedan f (0) = 12, s˚a kan kan vi best¨amma denna tidpunkt ur ekvationen f (t) = 6, dvs den ¨ar ett nollst¨alle f¨or funktionen 10e^−3t+2e^−5t−6 (t.ex. t ≈ 0.211327).

Vi skall d¨arf¨or studera iterativa metoder att best¨amma r¨otter, och b¨orjar med kvadratroten ur ett positivt tal a. Geometriskt kan problemet tolkas s˚a, att det g¨aller att konstruera en kvadrat, vars yta ¨ar a. Vi kan g¨ora det s˚a, att vi f¨orst gissar en rot, som vi t.ex. kallar x₀. En rektangel med sidorna x₀ och a/x₀ har d˚a ytan a. F¨or att f˚a rektangeln att mera likna en kvadrat, kan vi t.ex. r¨akna ut medeltalet av x0 och a/x₀: x₁ = ¹₂(x₀ + ^a

x0), och v¨alja detta till sida i en ny rektangel, och upprepa konstruktionen. Till slut f¨orvandlas rektangeln till en kvadrat (hoppas vi).

Vi kan l¨att skriva ett litet MATLAB–program som testar metoden:

>> a=5000;

>> x0=60;

>> for i=1:5

x1 = (x0 + a/x0)/2;

disp(sprintf(’%2.0f %17.14f’,i,x1)) x0 = x1;

end

1 71.66666666666666 2 70.71705426356590 3 70.71067840610468 4 70.71067811865476 5 70.71067811865476

Som synes, konvergerar den r¨att snabbt. Hur kan man konstruera ett bra utg˚angsv¨arde f¨or rotiterationerna?

Det ¨ar l¨att att se, att iterationsintervallet kan betydligt f¨orminskas. Om vi uttrycker a i formen a = p · 4^k,

1

4 < p < 1, d¨ar k ¨ar ett heltal, s˚a kan kvadratroten skrivas √

a = √

p · 2^k. Vi har allts˚a reducerat problemet till att ber¨akna kvadratroten f¨or ett tal inom intervallet [0.25, 1]. Ett bra utg˚angsv¨arde f¨or det reducerade kvadratrotsproblemet ¨ar d¨arf¨or L(p) = (1 + 2p)/3, eftersom denna funktion interpolerar f (p) = √

p i punkterna p = 0.25 och p = 1. Dessutom kan man visa, att |L(p) − √

p| ≤ 0.05 f¨or alla v¨arden av p inom detta intervall.

En uppskattning av felet efter k iterationer f˚ar vi ur ekvationen

x₁ − √

p = ¹₂



x₀ + p x₀



− √

p = ¹₂

x₀ − √

√ p x₀

² .

Antag, att x₀ = L(p) ¨ar utg˚angsv¨ardet och xk v¨ardet efter den k:te iterationen. Om felet ¨ar ek = x_k −√

p, s˚a f¨oljer av ekvationen ovan att e_k+1 = e²_k/(2x_k). Man kan ocks˚a visa, att uppskattningarna x_k alltid befinner sig inom intervallet [0.5, 1], s˚a att e_k+1 ≤ e²_k. H¨arav f¨oljer, att e₄ ≤ e²₃ ≤ e⁴₂ ≤ e⁸₁ ≤ e¹⁶₀ ≤ 0.05¹⁶, varf¨or fyra steg r¨acker till f¨or 16 siffors noggrannhet.

Programmet ovan kan d˚a modifieras p˚a f¨oljande s¨att (dock inte s˚a noggrannt f¨or sm˚a v¨arden av a)

>> a=5000;p=a; k=0;

>> while p>1 p=p/4;

k=k+1;

end

>> x0=(1+2*p)/3;

>> for i=1:5

x1=(x0+p/x0)/2; x=x1*2^k;

disp(sprintf(’%2.0f %17.14f’,i,x)) x0=x1;

end

1 70.73985496260360 2 70.71068413568874 3 70.71067811865501 4 70.71067811865476 5 70.71067811865476

En vanlig metod f¨or att s¨oka r¨otter ¨ar bisektionsmetoden. Den bygger p˚a satsen, att om en kontinuerlig funktion byter sitt f¨ortecken inom ett intervall, s˚a m˚aste den ha minst en rot inom intervallet. Detta resultat kan anv¨andas f¨or att begr¨ansa roten inom allt sn¨avare gr¨anser. Antag, att funktionens v¨arden i intervallets

¨andpunkter uppfyller villkoret f (a)f (b) ≤ 0 och att intervallets mittpunkt ¨ar m = (a + b)/2. D˚a g¨aller antingen f (a)f (m) ≤ 0 eller f (m)f (b) ≤ 0. I det f¨orstn¨amnda fallet vet vi att det finns en rot i intervallet [a, m], i det senare fallet vet vi att det finns en rot i intervallet [m, b]. Halveringsprocessen kan forts¨attas, tills vi n˚att en p˚a f¨orhand best¨amd toleransgr¨ans delta:

while abs(a-b) > delta if f(a)*f((a+b)/2) <= 0

b = (a+b)/2;

else

a = (a+b)/2;

end end

rot = (a+b)/2;

Programmet ¨ar n˚agot bristf¨alligt, f¨or det kan h¨anda, att while–slingan aldrig avslutas, om delta ¨ar mindre

¨an det tal, som anger r¨aknenoggrannheten. Vi kan korrigera detta genom att f¨or¨andra while–satsen till while abs(a-b) > delta + eps*max(abs(a), abs(b)).

D¨arigenom garanteras, att slingan avslutas, ¨aven om delta ¨ar f¨or litet. Dessutom anv¨ander programmet tv˚a funktionsber¨akningar vid varje iteration. F¨oljande program ¨ar en f¨orb¨attrad version av bisektionsmetoden:

function rot = bisect(fname,a,b,delta)

% Invariabler:

% fname funktionens (f(x)) namn

% a,b intervallet [a,b] d¨ar

% f ¨ar kontinuerlig, f(a)f(b) < 0

% delta icke-negativt reellt tal.

% Utvariabel:

% rot mittpunkten av ett intervall [a1,b1]

% f¨or vilket f(a1)f(b1)<=0 och

% |b1-a1| <= delta + eps*max(|a1|,|b1|) fa = feval(fname,a);

fb = feval(fname,b);

if fa*fb > 0

disp(’Roten ¨ar utanf¨or intervallet’) return

end

if nargin==3 delta = 0;

end

while abs(a-b) > delta+eps*max(abs(a),abs(b)) mid = (a+b)/2;

fmid = feval(fname,mid);

if fa*fmid<=0

% Det finns en rot innanf¨or [a,mid].

b = mid;

fb = fmid;

else

% Det finns en rot innanf¨or [mid,b].

a = mid;

fa = fmid;

end end

rot = (a+b)/2;

St¨orsta delen av tiden i programmet bisect ˚atg˚ar till att ber¨akna funktionsv¨arden.

Felet i bisektionsmetoden minskar med h¨alften vid varje steg. Om det k:te intervallet betecknas [ak, b_k], s˚a ¨ar |ak − b_k| ≤ |a₀ − b₀|/2^k varf¨or iterationsprocessen alltid konvergerar.

Om x_k = (a_k + b_k)/2 ¨ar den k:te approximationen till roten, s˚a g¨aller f¨or roten x_∗ villkoret

|xk − x∗| ≤ |ak − bk|

2 ≤ |a0 − b0| 2^k+1 .

Man s¨ager att en r¨acka xk konvergerar linj¨art mot x∗ om det finns en s˚adan konstant c, 0 ≤ c < 1 och ett heltal k₀ att |x_k+1 − x_∗| ≤ c|x_k − x_∗| f¨or alla k ≥ k₀ (c = 1/2 i detta fall). Vi kan till¨ampa bisect p˚a funktionen f (x) = tan(x/4) − 1 och s¨atta [a0, b₀] = [2, 4]:

a a_k b_k a_k - b_k

0 2.00000000000000 4.00000000000000 2.00000000000000 1 3.00000000000000 4.00000000000000 1.00000000000000 2 3.00000000000000 3.50000000000000 0.50000000000000 3 3.00000000000000 3.25000000000000 0.25000000000000 4 3.12500000000000 3.25000000000000 0.12500000000000 5 3.12500000000000 3.18750000000000 0.06250000000000 6 3.12500000000000 3.15625000000000 0.03125000000000

. . . .

43 3.14159265358967 3.14159265358990 0.00000000000023 44 3.14159265358978 3.14159265358990 0.00000000000011 45 3.14159265358978 3.14159265358984 0.00000000000006 46 3.14159265358978 3.14159265358981 0.00000000000003 47 3.14159265358978 3.14159265358980 0.00000000000001 48 3.14159265358979 3.14159265358980 0.00000000000001 49 3.14159265358979 3.14159265358980 0.00000000000000

Som vi ser, ¨ar konvergensen inte s¨arskilt snabb. Det beh¨ovs ca tre iterationer f¨or att ber¨akna en ny siffra i π.

Med Newtons metod kan man ber¨akna nollst¨allen betydligt snabbare. Antag, att vi k¨anner v¨ardet av en funktion f (x) och dess derivata f⁰(x) i en punkt x = x_c. Tangenten till kurvan i denna punkt

L(x) = f (x_c) + (x − x_c)f⁰(x_c)

kan uppfattas som en linj¨ar approximation f¨or kurvan i denna punkt (se figuren).

Nollst¨allet x0 f¨or L(x) kan d˚a ber¨aknas ur ekvationen

x₀ = x_c − f (x_c) f⁰(x_c).

Genom att upprepa denna formel f˚ar man en algoritm, som beskrivs genom f¨oljande MATLAB–program:

xc = input(’Ange begynnelsev¨ardet:’);

fc = feval(fname,xc);

dfc = feval(dfname,xc);

while input(’Nytt Newton-steg? (0=nej, 1=ja)’);

xnew = xc - fc/dfc;

xc = xnew;

fc = feval(fname,xc);

dfc = feval(dfname,xc);

end

Programmet f¨oruts¨atter, att fname och dfname ¨ar str¨angar, som inneh˚aller funktionsnamnet, resp.

derivatans namn. Om vi anv¨ander Newtons metod f¨or att ber¨akna nollst¨allet f¨or f (x) = tan(x/4) − 1 f˚ar vi en mycket snabbare konvergens ¨an med bisektionsmetoden:

k x_k |x_k - pi|

0 1.00000000000000 2.14159265358979 1 3.79631404657234 0.65472139298255 2 3.25943543617547 0.11784278258568 3 3.14513155420752 0.00353890061772 4 3.14159578639006 0.00000313280027 5 3.14159265359225 0.00000000000245 6 3.14159265358979 0.00000000000000

Som vi ser, verkar felet att bli kvadrerat efter ett visst antal iterationer. En s˚adan metod s¨ags konvergera kvadratiskt. I detta fall finns det ett heltal k₀ och en positiv konstant c som uppfyller villkoret

|x_k+1 − x_∗| ≤ c|x_k − x_∗|²

f¨or alla k ≥ k₀.

Newtons metod ¨ar dock inte alltid s¨arskilt stabil. Om t.ex. derivatans v¨arde i nollst¨allet f⁰(x_∗) ¨ar mycket litet, s˚a ¨ar tangenten n¨astan parallell med x–axeln, och det blir sv˚art att ber¨akna nollst¨allet noggrant. Om f⁰(x_c) ¨ar litet i f¨orh˚allande till f (x_c), s˚a kan korrektionen till roten bli alltf¨or stor, och Newton–steget f¨or oss alltf¨or l˚angt bort fr˚an nollst¨allet. Ett typiskt exempel ¨ar f (x) = arctan(x). Newton–korrektionen till roten ¨ar i detta fall x0 = x_c − (1 + x²_c) arctan(x_c). Man kan visa, att om |x_c| > 1.3917, s˚a ¨ar |x0| > |xc|. Detta betyder, att iterationerna divergerar om utg˚angsv¨ardet ¨ar utanf¨or intervallet [−1.3917, 1.3917]. Vi kan h¨arav dra den slutsatsen, att funktionen f inte kan vara alltf¨or olinj¨ar, och f⁰ inte alltf¨or n¨ara noll, om Newtons metod skall fungera.

Man kan visa, att om f⁰ inte ¨andrar f¨ortecken i n¨arheten av x∗, om f inte ¨ar alltf¨or olinj¨ar, och om Newton–

iterationerna p˚ab¨orjas tillr¨ackligt n¨ara ett nollst¨alle, s˚a garanteras en kvadratisk konvergens. Att avsluta en Newton–process ¨ar inte alldeles l¨att, eftersom man inte vet hur n¨ara minimet man kommit. En m¨ojlighet

¨ar att sluta d˚a |xk+1 − xk| ¨ar tillr¨ackligt litet, eftersom av lim xk = x_∗ f¨oljer lim |x_k+1 − xk| = 0.

Detta g¨aller dock inte omv¨ant. Om t.ex. f (x) = tan x och xc = π/2 − 10⁻⁵, s˚a ¨ar |x0 − xc| = tan x/(1 + tan² x) ≈ 10⁻⁵ fast¨an den n¨armaste roten ¨ar x = 0.

F¨or att bem¨astra dessa problem kan man kombinera Newtons metod med bisektionsmetoden. I b¨orjan av varje steg kontrolleras f¨orst att roten befinner sig innanf¨or intervallet [a, b], och att x_c ¨ar en av intervallets

¨

andpunkter. Om punkten

x₀ = x_c − f (x_c) f⁰(x_c)

h¨or till intervallet [a, b], s˚a ¨ar allt ok, och vi forts¨atter med att behandla antingen [a, x0] eller [x₀, b].

D¨arp˚a s¨atts xc lika med x₀. Om Newton–steget f¨or oss bort fr˚an intervallet [a, b], s˚a g¨or vi ett bisektions- steg, och s¨atter xc till (a + b)/2. F¨or att kontrollera om Newton–steget h˚alls innanf¨or intervallet [a, b]

anv¨ander vi funktionen

function ok = stegin(x,fx,dfx,a,b)

% Invariabler:

% x v¨ardet av x.

% fx v¨ardet av f i x.

% dfx v¨ardet av f’ i x.

% a,b anger intervallet [a,b]

% Utvariabel:

% ok 1 om Newton-steget x - fx/dfx h˚alls inom [a,b]

% 0 om inte.

if dfx > 0

ok = ((a-x)*dfx <= -fx) & (-fx <= (b-x)*dfx);

elseif dfx < 0

ok = ((a-x)*dfx >= -fx) & (-fx >= (b-x)*dfx);

else

ok = 0;

end

F¨or att f¨ors¨akra oss om att iterationerna tar slut, anv¨ands de tre f¨oljande konvergensvillkoren:

- L¨angden av intervallet som granskas skall vara mindre ¨an tolx, som ¨ar en best¨amd tolerans. Roten kommer d˚a inte att skilja sig fr˚an en riktig rot mer ¨an tolx.

- Absoluta v¨ardet av f (xc) ¨ar mindre eller lika med tolf, vilket inte beh¨over betyda att x_c ¨ar n¨ara en riktig rot.

- Antalet funktionsber¨akningar ¨overskrider ett positivt tal nEvalsMax. Detta inneb¨ar att de b˚ada tidigare n¨amnda villkoren inte ¨ar uppfyllda.

H¨ar ¨ar det slutliga programmet:

function [x,fx,nEvals,aF,bF] = newton(fName,dfName,a,b,tolx,tolf,nEvalsMax)

% Invariabler:

% fName namnet p˚a funktionen f(x).

% dfName namnet p˚a derivatafunktionen f’(x).

% a,b roten till f(x) s¨oks inom intervallet [a,b]

% och f(a)*f(b)<=0.

% tolx,tolf avslutningskriterier.

% nEvalsMax maximiantalet derivataber¨akningar.

% Utvariabler:

% x Ett approximativt nollst¨alle f¨or f.

% fx V¨ardet av f i x.

% nEvals Antalet derivataber¨akningar som beh¨ovdes.

% aF,bF Det slutliga intervallet [aF,bF].

fa = feval(fName,a);

fb = feval(fName,b);

if fa*fb>0

disp(’Roten inte innanf¨or intervallet’) return

end x = a;

fx = feval(fName,x);

dfx = feval(dfName,x);

disp(sprintf(’%20.15f %20.15f %20.15f’,a,x,b)) nEvals = 1;

while (abs(a-b)>tolx)&(abs(fx)>tolf)&((nEvals<nEvalsMax)|(nEvals==1))

%roten inom [a,b] och x = a eller x = b.

if stegin(x,fx,dfx,a,b)

%Ett Newton-steg:

disp(’Newton’)

x = x-fx/dfx;

else

%Ett bisektionssteg:

disp(’Bisektion’) x = (a+b)/2;

end

fx = feval(fName,x);

dfx = feval(dfName,x);

nEvals = nEvals+1;

if fa*fx<=0

% En rot inom [a,x]. V¨alj x till h¨oger ¨andpunkt.

b = x;

fb = fx;

else

% En rot inom [x,b]. V¨alj x till v¨anster ¨andpunkt.

a = x;

fa = fx;

end

disp(sprintf(’%20.15f %20.15f %20.15f’,a,x,b)) end

aF = a;

bF = b;

Detta program utf¨or normalt ett antal bisektionssteg innan Newton–iterationerna b¨orjar. H¨ar ¨ar resultatet, d˚a man ber¨aknar ett nollst¨alle f¨or f (x) = sin(x) inom intervallet [−7π/2, 15π + 0.1]:

Stegtyp a x b

-10.995574287564276 -10.995574287564276 47.223889803846896 Bisektion -10.995574287564276 18.114157758141312 18.114157758141312 Bisektion -10.995574287564276 3.559291735288517 3.559291735288517 Newton 3.115476144648328 3.115476144648328 3.559291735288517 Newton 3.115476144648328 3.141598592990409 3.141598592990409 Newton 3.141592653589793 3.141592653589793 3.141598592990409

Vi skall nu i korthet diskutera minimering av funktioner. Som ett exempel kan vi v¨alja v¨axelverkningspo- tentialen mellan Na⁺ och Cl⁻ jonerna i en NaCl–molekyl, som kan beskrivas av modellfunktionen

V (r) = − e²

4π₀r + αe^−r/ρ,

d¨ar e ¨ar elektronladdningen, och α = 1.09 · 10³eV och ρ = 0.330˚A ¨ar tv˚a parametrar, som beskriver v¨axelverkan mellan jonerna. Om vi substituerar x = r/ρ och s¨atter in talv¨arden, s˚a kan funktionen skrivas sortl¨ost i formen

f (x) = −0.04

x + e^−x ≡ V (r) α .

Denna funktion har ett minimum, som anger bindningsl¨angden f¨or molekylen NaCl. Detta ¨ar ett exempel p˚a ett endimensionellt optimeringsproblem. Funktionen f , som vi f¨ors¨oker optimera kallas objektivfunktionen.

Ett enkelt s¨att att studera en dylik funktion, ¨ar att upprita den och studera grafen f¨or att finna minimet.

Newtons metod, till¨ampad p˚a ekvationen f⁰(x) = 0, ¨ar en effektiv metod, ifall vi l¨att kan ber¨akna funktionens derivator, och vi k¨anner en god approximation f¨or minimet. Det finns ocks˚a enkla metoder att ber¨akna ett minimum f¨or en endimensionell funktion, d¨ar man inte beh¨over ber¨akna derivator.

En av dessa baserar sig p˚a det gyllene snittet och kan till¨ampas p˚a unimodala funktioner. En funktion f (x) s¨ags vara unimodal inom ett givet intervall [a, b], om det finns en punkt x_∗ inom detta intervall, s˚a att funktionen ¨ar str¨angt monotont avtagande inom [a, x∗] och str¨angt monotont v¨axande inom [x∗, b], dvs den har endast ett minimum inom intervallet.

Vi antar vidare, att vi har ber¨aknat funktionsv¨ardena i fyra punkter a1, a₂, a₃ och a₄, och att vi p˚a grund av detta vet att minimet befinner sig inom intervallet (a₁, a₄) (vi antar dessutom, att a₁ < a₂ < a₃ < a₄ g¨aller). Utg˚angsv¨ardena kan best¨ammas genom att man ber¨aknar funktionsv¨ardet i en viss punkt, och d¨arp˚a ger ett tillskott till argumentet tills funktionen b¨orjar v¨axa. Om detta inte sker, g˚ar man i motsatt riktning.

Antag nu, att dessa fyra a–v¨arden uppfyller ekvationerna a3 − a₁ = a₄ − a₂ = γ(a₄ − a₁), d¨ar γ = 2/(1 + √

5) ≈ 0.618034 . . .. Genom att testa funktionsv¨ardena f (a₁), f (a₂), . . . , f (a₄) ¨ar det m¨ojligt att f˚a reda p˚a inom vilket av de tv˚a lika stora intervallen (a₁, a₃) eller (a₂, a₄) minimet ligger.

Vi kan f¨or enkelhetens skull anta att det ligger inom intervallet (a₁, a₃). Funktionen ber¨aknas d¨arp˚a i en ny punkt a₅, som uppfyller villkoret a₅ − a1 = a₃ − a2. I detta fall g¨aller a3 − a5 = γ(a₃ − a1). Vi inser detta, om γ⁻¹ = γ + 1 till¨ampas i formeln ovan, varp˚a vi f˚ar a₄ − a1 = γ(a₃ − a1) + a₃ − a1, dvs γ(a₃ − a₁) = a₄ − a₃ = a₂ − a₁ = a₃ − a₅. Om vi nu betecknar punkterna a⁰₁ = a₁, a⁰₂ = a₅, a⁰₃ = a₂ och a⁰₄ = a₃, s˚a finner vi situationen motsvara utg˚angspunkten med undantag av att intervall¨angden reducerats med beloppet γ: a⁰₄ − a⁰₁ = γ(a₄ − a₁). Som av bilden framg˚ar, kommer allts˚a minimet att inneslutas mellan allt tr˚angare gr¨anser.

Gyllene snitt–metoden konvergerar linj¨art, liksom bisektionsmetoden. Newtons metod ¨ar snabb, men den till¨ampas p˚a derivatans nollst¨alle s˚a man beh¨over b˚ade f¨orsta och andra derivatan. MATLAB–funktionen fminbnd kombinerar gyllene snitt–metoden med en parabolisk metod, som inte beh¨over derivator.

Vi kan prova funktionen fminbnd p˚a potentialfunktionen f¨or NaCl molekylen:

function y = nacl(x)

y = -0.04./x + exp(-x);

Funktionen faller ganska brant, har ett maximum n¨ara 0.2, passerar x–axeln, och n¨armar sig x–axeln p˚a nytt fr˚an negativa sidan. Om minimet antas vara i intervallet [4, 10] n˚ar vi det efter 13 iterationer:

>> xmin=fminbnd(’nacl’,4,10,optimset(’TolX’,1e-14,’Display’,’iter’))

Func-count x f(x) Procedure

1 6.2918 -0.00450605 initial

2 7.7082 -0.00474015 golden

3 8.58359 -0.0044729 golden

4 7.40247 -0.00479386 parabolic

5 7.2692 -0.004806 parabolic

6 6.89587 -0.00478862 golden

7 7.1682 -0.00480949 parabolic

8 7.14704 -0.00480953 parabolic

9 7.15417 -0.00480955 parabolic

...

13 7.1543 -0.00480955 parabolic

Bindningsl¨angden f¨or NaCl ¨ar allts˚a r = 0.33 · 7.1543˚A = 2.3609˚A.