%% Rectangle inscribed inside an ellipse
% Find a rectangle inscribed inside an ellipse that maximizes 1) the circumference, 2) the area.
% The ellipse is given by the formula
% \frac{x_1^2}{a^2} + \frac{x_2^2}{b^2} = 1
%% Maximum circumference

syms x1 x2 a b lambda

f1 = 4*(x1+x2)
f1 = -f1

h = x1^2/a^2 + x2^2/b^2 - 1

L1 = f1 + lambda*h

gL1 = gradient(L1,[x1,x2,lambda])

%
%% Solving for the stationary point

sol1 = solve(gL1,[x1,x2,lambda])

x1opt = sol1.x1
x2opt = sol1.x2
lambdaopt = sol1.lambda

%
%% Plotting

A = 3; 
B = 1;

y1 = double(subs(x1opt(1),[a,b],[A,B]))
y2 = double(subs(x2opt(1),[a,b],[A,B]))

fimplicit(@(x1,x2) x1.^2/A^2+x2.^2/B^2-1); axis equal; hold on
hold on
rectangle('Position',[-y1 -y2 2*y1 2*y2])
hold off

%
%% Maximum area

f2 = 4*x1*x2
f2 = -f2

L2 = f2 + lambda*h

gL2 = gradient(L2,[x1,x2,lambda])

%
%% Solving for the stationary point

sol2 = solve(gL2,[x1,x2,lambda])

x1opt = sol2.x1
x2opt = sol2.x2
lambdaopt = sol2.lambda

%
%% Plotting

y1 = abs(double(subs(x1opt(1),[a,b],[A,B])))
y2 = abs(double(subs(x2opt(1),[a,b],[A,B])))

fimplicit(@(x1,x2) x1.^2/A^2+x2.^2/B^2-1); axis equal; hold on
hold on
rectangle('Position',[-y1 -y2 2*y1 2*y2])
hold off