% Script 2 for Week 3

The idea of the methods of finding a max or min of a continuous function f on an interval as outlined in Chapter 7 all center around picking 4 consecutive points x1-x2-x3-x4 for which your exploration of the function shows that there is definitely exactly ONE local extreme value (max OR min) inside (x1,x4). Your job is to localize that extreme value efficiently.

You start out with the brackets x1 and x4. The methods presented differ only in how x2 and x3 are obtained.

Suppose your algorithm is designed to find a local minimum of f.

The actual calculation inside each algorithm is to calculate f(x1), f(x2), f(x3) and f(x4). One of these four numbers is smallest. That smallest one is on the left half of the list f(x1)-f(x2)-f(x3)-f(x4) or on the right half.

If it is on the left half you retain x1-x2-x3 and use your method to pick a new number in (x1,x3) and now you are back to four numbers selected in a SMALLER interval [x1,x3] which I will relabel as newx1-newx2-newx3-newx4. 

If it is on the right half you retain x2-x3-x4 and use your method to pick a new number in (x2,x4) and now you are back to four numbers selected in a SMALLER interval [x2,x4] which I could relabel as newx1-newx2-newx3-newx4.

I have already calculated the function values at three of these points, and I will recycle these function calls. I just need to calculate one new function value.

Repeat till the interval is small enough for you, and the claim is that each step these smaller-and-smaller intervals capture the local extreme value you want.

But why is that?  There are 8 possible combinations of ``up-or-down'' on these three intervals.

We assume that other investigation has found that there IS a local minimum and No local maximum on the interval [x1,x4]. So there cannot be two local minima on [x1,x4].

1.  \ \ \   min in right segment (else max in first)
2.  / \ \   impossible: has a max inside
3.  \ / \   impossible: has a max inside
4.  \ \ /    min in right segment (else max in first)
5.  \ / /    min in left segment  (else max in second)
6.  / \ /   impossible: has a max inside
7.  / / \   impossible: has a max inside
8.  / / /    min in left segment  (else max in second)

So you proceed to
1. Throw away x1
4. Throw away x1
5. Throw away x4
8. Throw away x4
and the two segments that remain will always capture the minumum.


%**********************************************************

% IF YOU WANT AN EXTREME VALUE OR ROOT OF A FUNCTION, GRAPH IT FIRST!!!
 
 
 fxy=@(x) 4*x(1)+2*x(2)+x(1)^2-2*x(1)^4+2*x(1)*x(2)-3*x(2)^2;
 
[ x,y] = meshgrid([-3:.05:3]);
for k=1:length(x);
for w=1:length(y);
Z1(k,w) = fxy([x(k,w), y(k,w)]);
end
end

surfc(x,y,Z1)
colorbar
shading('interp')

fminsearch(fxy,[.85,.4])

g=@(x) -fxy(x)

[a b]=fminsearch(g,[.85,.4])

%**********************************************************


% Now some fun with matrices....

% We want to solve matrix equation Ax=B for variable "stack" x:

A=[2 1 0; 6 2 4; 0 2 4]
B=[4;7;1]

A\B

A^(-1)*B

G=[A B]

rref(G)

Tridiag([0 6 2],[2 2 4],[1 4 0],B)

tic
for k=1:10000
A\B;
end
toc

tic
for k=1:10000
A^(-1)*B;
end
toc

tic
for k=1:10000
rref(G);
end
toc

tic
for k=1:10000
Tridiag([0 6 2],[2 2 4],[1 4 0],B);
end
toc




% The script between the two lines below took about 45 minutes to run on my late-2013 Macbook Pro. If it takes too long on your computer you can change lim to a smaller interval, maybe lim=200:700. It plays a sound when it's done.

% Figure 6 displays the absolute time to solve a random matrix equation as a function of the matrix size using the backslash operator. Figure 7 suggests the squared dependency of this time. Figure 8 displays the ratio of time taken using the inverse versus backslash and figure 9 does the same for the rref solution. These last two suggest a size-cubed dependence of solution time on matrix dimension and that rref is definitely NOT optimized in MATLAB.

%************************************************************
A = cell(1,1000);
b = cell(1,1000);
for n = 1:1000;
    A{n} = rand(n,n);
    b{n} = rand(n,1);
    end
    
t = ones(1,1000);
s = ones(1,1000);
w = ones(1,1000);

lim=200:1000;

for n = lim
    tic;
    x = A{n}\b{n};
    t(n) = toc;
end

for n = lim
    tic;
    x = inv(A{n})*b{n};
    s(n) = toc;
end

for n = lim
    tic;
    x = rref([A{n} b{n}]);
    w(n) = toc;
end


figure(6)
plot(lim,t(lim))
title('Absolute Time, Backslash','FontSize',20,'FontWeight','bold')
ylabel('Time (Sec)')
xlabel('Matrix Size')

figure(7)
plot(lim, t(lim)./(lim.^2))
title('Ratio of Backslash Time to Squared Matrix Size','FontSize',20,'FontWeight','bold')
ylabel('Scaled Time')
xlabel('Matrix Size')

figure(8)
plot(lim, s(lim)./t(lim))
title('Ratio of Inverse to Backslash','FontSize',20,'FontWeight','bold')
ylabel('Relative Time')
xlabel('Matrix Size')

figure(9)
plot(lim, w(lim)./t(lim))
title('Ratio of RREF to Backslash','FontSize',20,'FontWeight','bold')
ylabel('Relative Time')
xlabel('Matrix Size')

load handel 
sound(y)

%****************************************