Application Center - Maplesoft

# Fibonacci search method for unimodal optimization

You can switch back to the summary page by clicking here.

fibsearch.mws

Fibonacci Search in Optimization of Unimodal Functions

by Dr. William P. Fox, Department of Mathematics
Francis Marion University, Florence, SC 29501, wfox@fmarion.edu
(c) 2002 William P. Fox

This program performs the Fibonacci Line Search algorithm to find the maximum of a unimodal function, f(x) , over an interval, a < x < b . The program calculates the number of iterations required to insure the final interval is within the user-specified tolerance. This is found by solving for the smallest value of n that makes this inequality true: Fn >(b-a)/tolerance, where n is the Fibonacci number from the sequence {F0,F1, F2, ...}. The Fibonacci search concept involves placing two experiments between [a,b] using the ratios of Fibonacci numbers. (The limit of the ratio of Fibonacci numbers is the golden section 0.618 but the Fibonacci method converges quicker.) One experiment is placed at position, , and the other at position, . The function, to be maximized, is evaluated at these two points and the functional values are compared. We want to keep the larger functional value (in our maximization problem) and its corresponding opposite end -point. At the end of the required iterations, the final interval is the answer. At times when the final answer must be a single point and not an interval, the convention of selecting the midpoint is provided. This program works when the function is not differentiable and we are looking for a solution. If you need to minimize a function, then multiple the function by (-1) and find the maximum.

To find a maximum solution to given a function, f(x), on the interval [a, b] where the function, f(x), is unimodal.

Overview of the algorithm

INPUT : endpoints a, b; tolerance, t, Fibonacci sequence
OUTPUT : final interval [ai, bi], f(midpoint)
Step 1 . Initialize the tolerance, t >0.
Step 2 . Set Fn > (b-a)/t as the smallest Fn and define the test points

Step 3
. Calculate f(x1) and f(x2)
Step 4 . Compare f(x1) and f(x2)
(a.) If f(x1) < f(x2), then the new interval is [x1, b]:
a becomes the previous x1
b does not change
x1 becomes the previous x2
n = n-1
Find the new x2 using the formula in Step 2.
(b.) If f(x1) > f(x2), then the new interval is [a, x2]:
a does not change
b becomes the previous x2
x2 becomes the previous x1
n = n-1
Find the new x1 using the formula in Step 2.
Step 5 . If the length of the new interval from Step 4 is less than the tolerance specified, the stop. Otherwise go back to Step 3.
Step 6 . Estimate x* as the midpoint of the final interval and compute, f(x*), the estimated maximum of the function.
STOP

To utilize the Fibonacci routine you need the following:

User inputs are:

The user enters the function f using

f:=x-><enter the expression in x>

Type FIBSearch(f, a, b, tolerance) for specific values of a, b, and the tolerance .

The output is the iterative process (each step).

The last output provided is the midpoint of the final interval and the value of f(x) at that point.

Program Code

> FIBSearch:=proc(f::procedure,a::numeric,b::numeric,T::numeric)

> local L, M, j, C, x1, x2, N, val;

> L:=(b-a)/T;

> M:=round(L);

> j:=0;

> label_2;

> C:=fib(j);

> if C<L then

> j:=j+1;

> goto(label_2);

> else

> N:=j;

> end if;

>

> x1:=evalf(a+(fib(N-2)/fib(N))*(b-a));

> x2:=evalf(a+(fib(N-1)/fib(N))*(b-a));

> printf("The interval [a,b] is [% 4.2f,% 4.2f]and user specified tolerance level is% 6.5f.\n",a,b,T);

> ### WARNING: %x or %X format should be %y or %Y if used with floating point arguments
printf("The first 2 experimental endpoints are x1= % 6.3f and x2 = % 6.3f. \n",x1,x2);

> printf(" \n");

> printf(" \n");

>

> printf(" Iteration x(1) x(2) f(x1) f(x2) Interval \n");

> iterate(f,a,b,N,x1,x2);

>

> val:=f(mdpt);

> printf(" \n");

> printf(" \n");

> printf("The midpoint of the final interval is% 9.6f and f(midpoint) = % 7.3f. \n",mdpt, val);

> printf(" \n");

> printf(" \n");

> ### WARNING: %x or %X format should be %y or %Y if used with floating point arguments
printf("The maximum of the function is % 7.3f and the x value = % 9.6f \n",fkeep,xkeep);

> printf(" \n");

> printf(" \n");

> end:

> iterate:=proc(f::procedure,a::numeric,b::numeric, N::posint,x1::numeric,x2::numeric)

> local x1n,x2n,an,bn,i,fx1,fx2,j,f1,f2,fmid;

> global mdpt,fkeep,xkeep,fib;

> i:=1;

> x1n(1):=x1;

> x2n(1):=x2;

> an(1):=a;

> bn(1):=b;

> i:=1;

> for j from 1 to N do

> fx1(i):=evalf(f(x1n(i)));

> fx2(i):=evalf(f(x2n(i)));

> if fx1(i)<=fx2(i) then

> an(i+1):=x1n(i);

> bn(i+1):=bn(i);

> x1n(i+1):=x2n(i);

> x2n(i+1):=an(i+1)+(fib(N-i-1)/fib(N-i))*(bn(i+1)-an(i+1));

> else

> an(i+1):=an(i);

> bn(i+1):=x2n(i);

> x2n(i+1):=x1n(i);

> x1n(i+1):=an(i+1)+(fib(N-i-2)/fib(N-i))*(bn(i+1)-an(i+1));

> fi;

> i:=i+1;

>

> printf(" % 3.0f % 11.4f % 10.4f % 10.4f %10.4f [% 6.4f, % 6.4f]\n",i,x1n(i),x2n(i),fx1(i-1),fx2(i-1),an(i),bn(i));

> mdpt := (an(i) + bn(i))/2;

>

>

> if ((i+2)=N) then

> if evalf(f(an(i))) > evalf(f(bn(i))) or evalf(f(an(i))) > evalf(f(mdpt)) then

>

> fkeep := f(an(i)); xkeep := an(i);

> else

> if evalf(f(bn(i))) > evalf(f(mdpt)) then

> fkeep := f(bn(i)); xkeep := bn(i);

> else

> fkeep := f(mdpt); xkeep := mdpt;

> end if;

> end if;

> end if;

> if(N-i-2)<0 then: goto(label_3):end if;

> end do;

> label_3;

> end:

>

> fib:=proc(n::numeric)

> option rememer;

> fib(0):=1:fib(1):=1;

> if n<2 then

> n;

> else

> fib(n-1)+fib(n-2);

> end if;

> end:

Examples

Example 1

Maximize the function f(x)=-x^2+21.6*x+3 over the interval [0,20] . . Via calculus the maximum is at 10.8. Let's see how the Fibonacci Search does.

> f:= x->-x^2+21.6*x+3;

> FIBSearch(f,0,20,.001);

`The interval [a,b] is [ 0.00, 20.00]and user specified tolerance level is .00100.`

`The first 2 experimental endpoints are x1=  7.639 and x2 =  12.361. `

` `

`  `

` Iteration    x(1)       x(2)      f(x1)      f(x2)           Interval `

`     2     12.3607    15.2786   109.6501   117.2043     [ 7.6393,  20.0000]`

`     3     10.5573    12.3607   117.2043    99.5818     [ 7.6393,  15.2786]`

`     4      9.4427    10.5573   119.5811   117.2043     [ 7.6393,  12.3607]`

`     5     10.5573    11.2461   117.7978   119.5811     [ 9.4427,  12.3607]`

`     6     10.1316    10.5573   119.5811   119.4410     [ 9.4427,  11.2461]`

`     7     10.5573    10.8204   119.1932   119.5811     [ 10.1316,  11.2461]`

`     8     10.8204    10.9830   119.5811   119.6396     [ 10.5573,  11.2461]`

`     9     10.7199    10.8204   119.6396   119.6065     [ 10.5573,  10.9830]`

`    10     10.8204    10.8825   119.6336   119.6396     [ 10.7199,  10.9830]`

`    11     10.7820    10.8204   119.6396   119.6332     [ 10.7199,  10.8825]`

`    12     10.7583    10.7820   119.6397   119.6396     [ 10.7199,  10.8204]`

`    13     10.7820    10.7967   119.6383   119.6397     [ 10.7583,  10.8204]`

`    14     10.7967    10.8057   119.6397   119.6400     [ 10.7820,  10.8204]`

`    15     10.7911    10.7967   119.6400   119.6400     [ 10.7820,  10.8057]`

`    16     10.7967    10.8002   119.6399   119.6400     [ 10.7911,  10.8057]`

`    17     10.8002    10.8022   119.6400   119.6400     [ 10.7967,  10.8057]`

`    18     10.7988    10.8002   119.6400   119.6400     [ 10.7967,  10.8022]`

`    19     10.8002    10.8009   119.6400   119.6400     [ 10.7988,  10.8022]`

`    20     10.7995    10.8002   119.6400   119.6400     [ 10.7988,  10.8009]`

`    21     10.8002    10.8002   119.6400   119.6400     [ 10.7995,  10.8009]`

`  `

`  `

`The midpoint of the final interval is 10.800154 and f(midpoint) =  119.640. `

`  `

`  `

`The maximum of the function is  119.640 and the x value =  10.799805 `

`  `

`  `

Example 2

Maximize a function whose derivative can not be solved explictly for 0.

Maximize the function, f(x) = x-exp(x) over the interval [-1,3].

> f:= x->x^2-exp(x);

> FIBSearch(f,0,3,.1);

`The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is .10000.`

`The first 2 experimental endpoints are x1=  1.147 and x2 =  1.853. `

` `

`  `

` Iteration    x(1)       x(2)      f(x1)      f(x2)           Interval `

`     2      1.8529     2.2941     1.1482     1.1937     [ 1.1471,  3.0000]`

`     3      2.2941     2.5588     1.1937     1.2027     [ 1.8529,  3.0000]`

`     4      2.5588     2.7353     1.2027     1.2036     [ 2.2941,  3.0000]`

`     5      2.4706     2.5588     1.2036     1.2028     [ 2.2941,  2.7353]`

`     6      2.3824     2.4706     1.2036     1.2036     [ 2.2941,  2.5588]`

`     7      2.4706     2.4706     1.2033     1.2036     [ 2.3824,  2.5588]`

`  `

`  `

`The midpoint of the final interval is 2.470588 and f(midpoint) =   1.204. `

`  `

`  `

`The maximum of the function is   1.204 and the x value =  2.558824 `

`  `

`  `

Example 3.

Maximizing a function that does not have a derivative.

> f:= x->-(abs(2-x)+abs(5-4*x)+abs(8-9*x));

> FIBSearch(f,0,3,.1);

`The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is .10000.`

`The first 2 experimental endpoints are x1=  1.147 and x2 =  1.853. `

` `

`  `

` Iteration    x(1)       x(2)      f(x1)      f(x2)           Interval `

`     2       .7059     1.1471    -3.5882   -11.2353     [ 0.0000,  1.8529]`

`     3      1.1471     1.4118    -5.1176    -3.5882     [ .7059,  1.8529]`

`     4       .9706     1.1471    -3.5882    -5.9412     [ .7059,  1.4118]`

`     5       .8824      .9706    -2.8824    -3.5882     [ .7059,  1.1471]`

`     6       .7941      .8824    -2.6471    -2.8824     [ .7059,  .9706]`

`     7       .8824      .8824    -3.8824    -2.6471     [ .7941,  .9706]`

`  `

`  `

`The midpoint of the final interval is  .882353 and f(midpoint) =  -2.647. `

`  `

`  `

`The maximum of the function is  -2.882 and the x value =   .970588 `

`  `

`  `

Example 5.

Minimize the function, . To minimize , we maximize . We keep the interval the same.

> f:= x->-2*x^2+4*x;

> FIBSearch(f,-1,2,.1);

`The interval [a,b] is [-1.00, 2.00]and user specified tolerance level is .10000.`

`The first 2 experimental endpoints are x1=   .147 and x2 =   .853. `

` `

`  `

` Iteration    x(1)       x(2)      f(x1)      f(x2)           Interval `

`     2       .8529     1.2941      .5450     1.9567     [ .1471,  2.0000]`

`     3       .5882      .8529     1.9567     1.8270     [ .1471,  1.2941]`

`     4       .8529     1.0294     1.6609     1.9567     [ .5882,  1.2941]`

`     5      1.0294     1.1176     1.9567     1.9983     [ .8529,  1.2941]`

`     6       .9412     1.0294     1.9983     1.9723     [ .8529,  1.1176]`

`     7      1.0294     1.0294     1.9931     1.9983     [ .9412,  1.1176]`

`  `

`  `

`The midpoint of the final interval is 1.029412 and f(midpoint) =   1.998. `

`  `

`  `

`The maximum of the function is   2.000 and the x value =   .985294 `

`  `

`  `

> f:=(x)->1-exp(-x)+1/(1+x);

> FIBSearch(f,0,20,.1);

`The interval [a,b] is [ 0.00, 20.00]and user specified tolerance level is .10000.`

`The first 2 experimental endpoints are x1=  7.639 and x2 =  12.361. `

` `

`  `

` Iteration    x(1)       x(2)      f(x1)      f(x2)           Interval `

`     2      4.7210     7.6395     1.1153     1.0748     [ 0.0000,  12.3605]`

`     3      2.9185     4.7210     1.1659     1.1153     [ 0.0000,  7.6395]`

`     4      1.8026     2.9185     1.2012     1.1659     [ 0.0000,  4.7210]`

`     5      2.9185     3.6052     1.1919     1.2012     [ 1.8026,  4.7210]`

`     6      2.4893     2.9185     1.2012     1.1900     [ 1.8026,  3.6052]`

`     7      2.2318     2.4893     1.2036     1.2012     [ 1.8026,  2.9185]`

`     8      2.4893     2.6609     1.2021     1.2036     [ 2.2318,  2.9185]`

`     9      2.4034     2.4893     1.2036     1.2033     [ 2.2318,  2.6609]`

`    10      2.4893     2.5751     1.2034     1.2036     [ 2.4034,  2.6609]`

`    11      2.4893     2.4893     1.2036     1.2036     [ 2.4034,  2.5751]`

`  `

`  `

`The midpoint of the final interval is 2.489270 and f(midpoint) =   1.204. `

`  `

`  `

`The maximum of the function is   1.203 and the x value =  2.403433 `

`  `

`  `