Technology for Chapter 7 and Chapter 13: Continuous Optimization (Numerical Techniques)
Chapter 7: Discrete Single Variable Techniques
Lab covers the following techniques:
Dichotomous
Golden Section
Maple:
Programs were written for each of these techniques.
DICHOTOMOUS SEARCH ALGORITHM- Maximization
Dr. William P. Fox, Naval PostgraduateSchool, Monterey, CA93953
> restart;
> DICHOTOMOUS:=proc(f::procedure,a::numeric,b::numeric,T::numeric,Ep::numeric)
> local x1,x2;
> x1:=(a+b)/2-Ep;
> x2:=(a+b)/2+Ep;
> 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");
> N:=ceil((ln(T/(b-a))/ln(0.5)));
> printf(" Iteration x(1) x(2) f(x1) f(x2) Interval \n");
> iterate(f,a,b,N,x1,x2,Ep);
> 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,Ep::numeric)
> local x1n,x2n,an,bn,i,fx1,fx2,j,f1,f2,fmid;
> global mdpt,fkeep,xkeep;
> i:=1;
> x1n(1):=x1;
> x2n(1):=x2;
> an(1):=a;
> bn(1):=b;
> i:=1;
> for j from 1 to N+1 do
> fx1(i):=f(x1n(i));
> fx2(i):=f(x2n(i));
> if fx1(i)<=fx2(i) then
> an(i+1):=x1n(i);
> bn(i+1):=bn(i);
> x1n(i+1):=(an(i+1)+bn(i+1))/2-Ep;
> x2n(i+1):=(an(i+1)+bn(i+1))/2+Ep;
> else
> an(i+1):=an(i);
> bn(i+1):=x2n(i);
> x1n(i+1):=(an(i+1)+bn(i+1))/2-Ep;
> x2n(i+1):=(an(i+1)+bn(i+1))/2+Ep;
> fi;
> printf(" % 3.0f % 11.4f % 10.4f % 10.4f %10.4f [% 6.4f, % 6.4f]\n",i,x1n(i),x2n(i),fx1(i),fx2(i),an(i),bn(i));
> mdpt := (an(i) + bn(i))/2;
> i:=i+1;
> if (i=N+1) then
> if (f(an(i)) > f(bn(i)) or f(an(i)) > f(mdpt)) then
> fkeep := f(an(i)); xkeep := an(i);
> else
> if (f(bn(i)) > f(mdpt)) then
> fkeep := f(bn(i)); xkeep := bn(i);
> else
> fkeep := f(mdpt); xkeep := mdpt;
> fi;
> fi;
> fi;
> od;
> end:
Warning, `N` is implicitly declared local to procedure `DICHOTOMOUS`
Warning, `val` is implicitly declared local to procedure `DICHOTOMOUS`
> f:=x->-x^2-2*x;
> DICHOTOMOUS(f,-3,6,.2,.01);
The interval [a,b] is [-3.00, 6.00]and user specified tolerance level is 0.20000.
The first 2 experimental endpoints are x1= 1.490 and x2 = 1.510.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 1.4900 1.5100 -5.2001 -5.3001 [-3.0000, 6.0000]
2 -0.7550 -0.7350 0.9400 0.9298 [-3.0000, 1.5100]
3 -1.8775 -1.8575 0.2300 0.2647 [-3.0000, -0.7350]
4 -1.3162 -1.2962 0.9000 0.9122 [-1.8775, -0.7350]
5 -1.0356 -1.0156 0.9987 0.9998 [-1.3162, -0.7350]
6 -0.8953 -0.8753 0.9890 0.9845 [-1.0356, -0.7350]
7 -0.9655 -0.9455 0.9988 0.9970 [-1.0356, -0.8753]
The midpoint of the final interval is-0.955469 and f(midpoint) = 0.998.
The maximum of the function is 0.999 and the x value = -1.035625
> f:=x->3*exp(.1*x)+cos(x)+sqrt(x);
> DICHOTOMOUS(f,0,2.5,.1,.01);
The interval [a,b] is [ 0.00, 2.50]and user specified tolerance level is .10000.
The first 2 experimental endpoints are x1= 1.240 and x2 = 1.260.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 1.2400 1.2600 4.8344 4.8312 [ 0.0000, 2.5000]
2 .6200 .6400 4.7932 4.8004 [ 0.0000, 1.2600]
3 .9300 .9500 4.8546 4.8553 [ .6200, 1.2600]
4 1.0850 1.1050 4.8524 4.8508 [ .9300, 1.2600]
5 1.0075 1.0275 4.8557 4.8553 [ .9300, 1.1050]
6 .9688 .9888 4.8557 4.8559 [ .9300, 1.0275]
The midpoint of the final interval is .998125 and f(midpoint) = 4.856.
The maximum of the function is 4.855 and the x value = .930000
> f:=x->-(x^2)-2*x;
> DICHOTOMOUS(f,-3,6,.2,.01);
The interval [a,b] is [-3.00, 6.00]and user specified tolerance level is .20000.
The first 2 experimental endpoints are x1= 1.490 and x2 = 1.510.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 1.4900 1.5100 -5.2001 -5.3001 [-3.0000, 6.0000]
2 -.7550 -.7350 .9400 .9298 [-3.0000, 1.5100]
3 -1.8775 -1.8575 .2300 .2647 [-3.0000, -.7350]
4 -1.3163 -1.2963 .9000 .9122 [-1.8775, -.7350]
5 -1.0356 -1.0156 .9987 .9998 [-1.3163, -.7350]
6 -.8953 -.8753 .9890 .9845 [-1.0356, -.7350]
7 -.9655 -.9455 .9988 .9970 [-1.0356, -.8753]
The midpoint of the final interval is -.990547 and f(midpoint) = 1.000.
The maximum of the function is .999 and the x value = -1.035625
Example 1
Maximize the function f(x)=1-exp(x)+1/(1+x) over the interval [0,20]..
> f:= x->1-exp(-x)+(1/(1+x));
> DICHOTOMOUS(f,0,20,.001,.001);
The interval [a,b] is [ 0.00, 20.00]and user specified tolerance level is .00100.
The first 2 experimental endpoints are x1= 9.999 and x2 = 10.001.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 9.9990 10.0010 1.0909 1.0909 [ 0.0000, 20.0000]
2 4.9995 5.0015 1.1599 1.1599 [ 0.0000, 10.0010]
3 2.4998 2.5018 1.2036 1.2036 [ 0.0000, 5.0015]
4 3.7496 3.7516 1.1870 1.1870 [ 2.4998, 5.0015]
5 3.1247 3.1267 1.1985 1.1985 [ 2.4998, 3.7516]
6 2.8122 2.8142 1.2022 1.2022 [ 2.4998, 3.1267]
7 2.6560 2.6580 1.2033 1.2033 [ 2.4998, 2.8142]
8 2.5779 2.5799 1.2036 1.2036 [ 2.4998, 2.6580]
9 2.5388 2.5408 1.2036 1.2036 [ 2.4998, 2.5799]
10 2.5193 2.5213 1.2036 1.2036 [ 2.4998, 2.5408]
11 2.5095 2.5115 1.2036 1.2036 [ 2.4998, 2.5213]
12 2.5144 2.5164 1.2036 1.2036 [ 2.5095, 2.5213]
13 2.5120 2.5140 1.2036 1.2036 [ 2.5095, 2.5164]
14 2.5107 2.5127 1.2036 1.2036 [ 2.5095, 2.5140]
15 2.5113 2.5133 1.2036 1.2036 [ 2.5107, 2.5140]
16 2.5117 2.5137 1.2036 1.2036 [ 2.5113, 2.5140]
The midpoint of the final interval is 2.512803 and f(midpoint) = 1.204.
The maximum of the function is 1.204 and the x value = 2.512346
Example 2.
>
> f:= x->-(x^2)-1;
> DICHOTOMOUS(f,-1,0.75,.25,.01);
The interval [a,b] is [-1.00, .75]and user specified tolerance level is .25000.
The first 2 experimental endpoints are x1= -.135 and x2 = -.115.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 -.1350 -.1150 -1.0182 -1.0132 [-1.0000, .7500]
2 .2975 .3175 -1.0885 -1.1008 [-.1350, .7500]
3 .0813 .1013 -1.0066 -1.0103 [-.1350, .3175]
4 -.0269 -.0069 -1.0007 -1.0000 [-.1350, .1013]
The midpoint of the final interval is .037188 and f(midpoint) = -1.001.
The maximum of the function is -1.018 and the x value = -.135000
Example 3.
> p := piecewise(x<=2,x/2,x>2,-x+3);
> f:= x->piecewise(x<=2,x/2,x>2,-x+3);
> plot(piecewise(x<=2,x/2,x>2,-x+3),x=0..3);
> DICHOTOMOUS(f,0,3,.25,.01);
1 1.4900 1.5100 The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is .25000.
The first 2 experimental endpoints are x1= 1.490 and x2 = 1.510.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 1.4900 1.5100 .7450 .7550 [ 0.0000, 3.0000]
2 2.2350 2.2550 .7650 .7450 [ 1.4900, 3.0000]
3 1.8625 1.8825 .9313 .9413 [ 1.4900, 2.2550]
4 2.0488 2.0688 .9513 .9313 [ 1.8625, 2.2550]
5 1.9556 1.9756 .9778 .9878 [ 1.8625, 2.0688]
The midpoint of the final interval is 2.012188 and f(midpoint) = .988.
The maximum of the function is .931 and the x value = 1.862500
Example 4.
Maximizing a function that does not have a derivative.
> f:= x->-(abs(2-x)+abs(5-4*x)+abs(8-9*x));
> DICHOTOMOUS(f,0,3,.1,.01);
The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is .10000.
The first 2 experimental endpoints are x1= 1.490 and x2 = 1.510.
Iteration x(1) x(2) f(x1) f(x2) Interval
1 1.4900 1.5100 .7450 .7550 [ 0.0000, 3.0000]
2 2.2350 2.2550 .7650 .7450 [ 1.4900, 3.0000]
3 1.8625 1.8825 .9313 .9413 [ 1.4900, 2.2550]
4 2.0488 2.0688 .9513 .9313 [ 1.8625, 2.2550]
5 1.9556 1.9756 .9778 .9878 [ 1.8625, 2.0688]
6 2.0022 2.0222 .9978 .9778 [ 1.9556, 2.0688]
The midpoint of the final interval is 1.988906 and f(midpoint) = .994.
The maximum of the function is .983 and the x value = 1.965625
Golden Section
GOLDEN SECTION SEARCH ALGORITHM-
This program performs the Golden Section 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 k that makes this inequality true: (b-a)(0.618^k)< tolerance. The Golden section concept involves placing two experiments between [a,b] using the Golden section ratios. One experiment is placed at position, a+0.382(b-a), and the other at position, a+0.618(b-a). 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 utilizie the Golden section routine you need the following:
User INPUTs-
The user enters the function f using
f:=x-<enter the expression in x>
Type GOLD(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.
Dr. William P. Fox, Department of Mathematics, FrancisMarionUniversity, Florence, SC29501
Dr. Margie Witherspoon, Department of Computer Science, FrancisMarionUniversity, Florence, SC29501
> restart;
> GOLD:=proc(f::procedure,a::numeric,b::numeric,T::numeric)
> local x1,x2;
> x1:=a+0.382*(b-a);
> x2:=a+0.618*(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");
> N:=ceil((ln(T/(b-a))/ln(0.618)));
> printf(" Iteration Interval x(1) x(2) f(x1) f(x2) \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;
> i:=1;
> x1n(1):=x1;
> x2n(1):=x2;
> an(1):=a;
> bn(1):=b;
> i:=1;
> for j from 1 to N+1 do
> fx1(i):=f(x1n(i));
> fx2(i):=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)+.618*(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)+.382*(bn(i+1)-an(i+1));
> fi;
> i:=i+1;
> printf(" % 3.0f [% 6.4f, % 6.4f] % 11.4f % 10.4f % 10.4f %10.4f \n",i-1, an(i-1),bn(i-1),x1n(i-1),x2n(i-1),fx1(i-1),fx2(i-1));
> mdpt := (an(i) + bn(i))/2;
> ##printf(" % 11.4f % 11.4f\n",an(i),bn(i));
> if (i=N) then
> if (f(an(i-1)) > f(bn(i-1)) or f(an(i-1)) > f(mdpt)) then
>
> fkeep := f(an(i)-1); xkeep := an(i-1);
> else
> if (f(bn(i-1)) > f(mdpt)) then
> fkeep := f(bn(i-1)); xkeep := bn(i-1);
> else
> fkeep := f(mdpt); xkeep := mdpt;
> fi;
> fi;
> fi;
> od;
> end:
Warning, `N` is implicitly declared local to procedure `GOLD`
Warning, `val` is implicitly declared local to procedure `GOLD`
The interval [a,b] is [ 0.00, 25.00]and user specified tolerance level is 0.25000.
The first 2 experimental endpoints are x1= 9.550 and x2 = 15.450.
Iteration Interval x(1) x(2) f(x1) f(x2)
1 [ 0.0000, 25.0000] 9.5500 15.4500 -66.3275 -381.3875
2 [ 0.0000, 15.4500] 5.9019 9.5500 23.9838 -66.3275
3 [ 0.0000, 9.5500] 3.6481 5.9019 39.8731 23.9838
4 [ 0.0000, 5.9019] 2.2545 3.6481 34.4491 39.8731
5 [ 2.2545, 5.9019] 3.6481 4.5086 39.8731 37.4033
6 [ 2.2545, 4.5086] 3.1156 3.6481 39.1760 39.8731
7 [ 3.1156, 4.5086] 3.6481 3.9765 39.8731 39.4548
8 [ 3.1156, 3.9765] 3.4444 3.6481 39.8074 39.8731
9 [ 3.4444, 3.9765] 3.6481 3.7732 39.8731 39.7900
10 [ 3.4444, 3.7732] 3.5700 3.6481 39.8773 39.8731
11 [ 3.4444, 3.6481] 3.5222 3.5700 39.8619 39.8773
The midpoint of the final interval is 3.585170 and f(midpoint) = 39.879.
The maximum of the function is 35.874 and the x value = 3.444442
Maximize f(x) = -|2-x|-|5-4x|-|8-9x| over the interval 0 < x < 3.
Example 3.
> f:= x->-3*x^2+21.6*x+1;
> GOLD(f,0,25,.1);
The interval [a,b] is [ 0.00, 25.00]and user specified tolerance level is 0.10000.
The first 2 experimental endpoints are x1= 9.550 and x2 = 15.450.
Iteration Interval x(1) x(2) f(x1) f(x2)
1 [ 0.0000, 25.0000] 9.5500 15.4500 -66.3275 -381.3875
2 [ 0.0000, 15.4500] 5.9019 9.5500 23.9838 -66.3275
3 [ 0.0000, 9.5500] 3.6481 5.9019 39.8731 23.9838
4 [ 0.0000, 5.9019] 2.2545 3.6481 34.4491 39.8731
5 [ 2.2545, 5.9019] 3.6481 4.5086 39.8731 37.4033
6 [ 2.2545, 4.5086] 3.1156 3.6481 39.1760 39.8731
7 [ 3.1156, 4.5086] 3.6481 3.9765 39.8731 39.4548
8 [ 3.1156, 3.9765] 3.4444 3.6481 39.8074 39.8731
9 [ 3.4444, 3.9765] 3.6481 3.7732 39.8731 39.7900
10 [ 3.4444, 3.7732] 3.5700 3.6481 39.8773 39.8731
11 [ 3.4444, 3.6481] 3.5222 3.5700 39.8619 39.8773
12 [ 3.5222, 3.6481] 3.5700 3.6000 39.8773 39.8800
13 [ 3.5700, 3.6481] 3.6000 3.6183 39.8800 39.8790
The midpoint of the final interval is 3.594161 and f(midpoint) = 39.880.
The maximum of the function is 39.879 and the x value = 3.585170
Example 2.
Maximizing a function that does not have a derivative.
> f:= x->-(abs(2-x)+abs(5-4*x)+abs(8-9*x));
> GOLD(f,0,3,.2);
The interval [a,b] is [ 0.00, 3.00]and user specified tolerance level is 0.20000.
The first 2 experimental endpoints are x1= 1.146 and x2 = 1.854.
Iteration Interval x(1) x(2) f(x1) f(x2)
1 [ 0.0000, 3.0000] 1.1460 1.8540 -3.5840 -11.2480
2 [ 0.0000, 1.8540] 0.7082 1.1460 -5.0848 -3.5840
3 [ 0.7082, 1.8540] 1.1460 1.4163 -3.5840 -5.9958
4 [ 0.7082, 1.4163] 0.9787 1.1460 -2.9149 -3.5840
5 [ 0.7082, 1.1460] 0.8755 0.9787 -2.7436 -2.9149
6 [ 0.7082, 0.9787] 0.8116 0.8755 -3.6382 -2.7436
7 [ 0.8116, 0.9787] 0.8755 0.9149 -2.7436 -2.6594
The midpoint of the final interval is 0.927087 and f(midpoint) = -2.708.
The maximum of the function is -3.191 and the x value = 0.843473
Example 3.
Minimize the function, 2x^2-4x. To minimize f(x), we maximize -f(x)= -2 x^2+4x. We keep the interval the same.
> f:= x->-2*x^2+4*x;
> GOLD(f,-1,2,.4);
The interval [a,b] is [-1.00, 2.00]and user specified tolerance level is 0.40000.
The first 2 experimental endpoints are x1= 0.146 and x2 = 0.854.
Iteration Interval x(1) x(2) f(x1) f(x2)
1 [-1.0000, 2.0000] 0.1460 0.8540 0.5414 1.9574
2 [ 0.1460, 2.0000] 0.8540 1.2918 1.9574 1.8297
3 [ 0.1460, 1.2918] 0.5837 0.8540 1.6534 1.9574
4 [ 0.5837, 1.2918] 0.8540 1.0213 1.9574 1.9991
5 [ 0.8540, 1.2918] 1.0213 1.1245 1.9991 1.9690
6 [ 0.8540, 1.1245] 0.9573 1.0213 1.9964 1.9991
The midpoint of the final interval is 1.040945 and f(midpoint) = 1.997.
The maximum of the function is 1.989 and the x value = 1.072886
>
Chapter 13
Steepest Ascent
NLP with Steepest Ascent
Bill Fox, Naval PostgraduateSchool
Computer Location
An Example of Optimization
The Problem
Consider a small company that is planning to install a central computer with cable links to five departments. According to their floor plan, the peripheral computers for the five departments will be situated as shown by the dark circles in Figure 1. The company wishes to locate the central computer so that the minimal amount of cable will be used to link to the five peripheral computers. Assuming that cable may be strung over the ceiling panels in a straight line from a point above any peripheral to a point above the central computer, the distance formula may be used to determine the length of cable needed to connect any peripheral to the central computer. Ignore all lengths of cable from the computer itself to a point above the ceiling panel immediately over that computer. That is, work only with lengths of cable strung over the ceiling panels.
The coordinates of the locations of the five peripheral computers are listed in Table 1.
X Y
15 60
25 90
60 75
75 60
80 25
Table 1. Grid Coordinates of Five Departments
> x:=[15,25,60,75,80];
> y:=[60,90,75,60,25];
> xy:=array(1..2,1..5,[x,y]);
> with(plots):
Warning, the name display has been redefined
> datalv:={seq([x[i],y[i]],i=1..5)};
> plot(datalv,labels=[X,Y],view=[0..100,0..100], style=point,symbol=box,color=black, thickness=3);
Solution and 3D graphical representation of function
> restart;
> with(plots):with(linalg):
Warning, the name changecoords has been redefined
Warning, the protected names norm and trace have been redefined and unprotected
> Dist(x,y):=(sqrt((x-15)^2+(y-60)^2)+sqrt((x-25)^2+(y-90)^2)+sqrt((x-60)^2+(y-75)^2)+sqrt((x-80)^2+(y-25)^2)+sqrt((y-60)^2+(x-75)^2));
> plot3d(Dist(x,y),x=0..100,y=0..100);
> distx:=diff(Dist(x,y),x);
> disty:=diff(Dist(x,y),y);
> with(plottools): with(plots):
> coolplot:= proc(k,y,x,a,b)
local tf;
> p:= plot3d(k,x=a..b,y=a..b);
> q:= contourplot(k,x=a..b,y=a..b,contours=50,axes=BOXED,color=BLACK);
> tf:= transform( (x,y) -> [x,y,-2.5] );
> display({p,tf(q)});
> end;
Warning, the name arrow has been redefined
Warning, the name arrow has been redefined
Warning, `p` is implicitly declared local to procedure `coolplot`
Warning, `q` is implicitly declared local to procedure `coolplot`
> k:=Dist(x,y);
> coolplot(k,y,x,50,80);
Gradient Search
Gradient Method-Steepest Ascent
Let's use a gradient search method to illustrate.
The algorithm inputs are:
n = maximum number of iterations
tol = tolerance to stop when magnitude of gradient is less than tolerance
ptx1, ptx2 are the initial (x0,y0) point
f is the multivariable function to be maximized
OUTPUTs are the approximate optimal point, (x,y), and the functional value at that point, f(x,y).
If you need to minimize a function, multiply f by -1 and use this algorithm.
Dr. William P. Fox, Department of Mathematics, FrancisMarionUniversity, Florence, SC29501 ()
Dr. William H. Richardson, Department of Mathematics, FrancisMarionUniversity, Florence, SC29501 ()
> restart:
> with(linalg):
Warning, the protected names norm and trace have been redefined and unprotected
> UP:=proc(n,tol,ptx1,ptx2,f)
> local numIter,numfeval,numgeval,f1,p1,p2,rv,x1pt,x2pt,temp,max,mg,v1,v2,newt,lam,nv1,nv2,Fvalue;
> f1:=f;temp:=grad(f1,vector([x1,x2]));
> p1 := unapply(temp[1],x1,x2);p2 := unapply(temp[2],x1,x2);
> x1pt:=ptx1;
> x2pt:=ptx2;
> rv:=vector([p1(x1pt,x2pt),p2(x1pt,x2pt)]):
> numIter:=1;numgeval:=1;numfeval:=1;
> printf("\n\n------");
> printf("------");
> printf("\n\n Initial Condition: ");
> printf("(%8.4f,%8.4f)\n\n",x1pt,x2pt);
> printf(" Iter Gradient Vector G magnitude G");
> printf(" x[k] Step Length\n\n");
> label_7;
> rv:=vector([p1(x1pt,x2pt),p2(x1pt,x2pt)]):
> numgeval:=numgeval+1;
> printf("%5d (%8.4f,%8.4f)",numIter,rv[1],rv[2]);
> max:=n;
> mg:=convert(sqrt(dotprod(rv,rv)),float);
> printf("%12.4f",mg);
> if(mg<tol or numIter>=max) then
> goto(label_6);
> else
> numIter:=numIter+1;
> fi;
> v1:=x1pt+t*rv[1];
> v2:=x2pt+t*rv[2];
> newt:=evalf(subs({x1=v1,x2=v2},f1));
> numfeval:=numfeval+1;
> lam:=fsolve(diff(newt,t)=0,t,maxsols=1);
> nv1:=evalf(subs({t=lam},v1));
> nv2:=evalf(subs({t=lam},v2));
> printf(" (%8.4f,%8.4f)%13.4f\n",x1pt,x2pt,lam);
> x1pt:=nv1;
> x2pt:=nv2;
> goto(label_7);
> label_6;
> printf("\n\n------");
> printf("------");
> printf("\n\n Approximate Solution: ");
> printf(" (%8.4f,%8.4f)\n",x1pt,x2pt);
> Fvalue:=evalf(subs(x1=x1pt,x2=x2pt,f));
> printf(" Maximum Functional Value: ");
> printf("%21.4f",Fvalue);
> printf("\n Number gradient evaluations:");
> printf("%22d",numgeval);
> printf("\n Number function evaluations:");
> printf("%22d",numfeval);
> printf("\n\n------");
> printf("------");
> end:
> f:=(sqrt((x1-15)^2+(x2-60)^2)+sqrt((x1-25)^2+(x2-90)^2)+sqrt((x1-60)^2+(x2-75)^2)+sqrt((x1-80)^2+(x2-25)^2)+sqrt((x2-60)^2+(x1-75)^2));
> UP(200,.005,0,3,f);
------
Initial Condition: ( 0.0000, 3.0000)
Iter Gradient Vector G magnitude G x[k] Step Length
1 ( -2.9312, -3.5666) 4.6166 ( 0.0000, 3.0000) -18.6620
2 ( -.3001, .2466) .3884 ( 54.7025, 69.5607) -6.6799
3 ( -.0133, -.0161) .0209 ( 56.7069, 67.9134) -9.2843
4 ( .0018, -.0015) .0023
------
Approximate Solution: ( 56.8299, 68.0631)
Maximum Functional Value: 157.6635
Number gradient evaluations: 5
Number function evaluations: 4
------