function [R,H] = sampleproblem5(area,volume,jacflag,R0,H0); %[R,H] = sampleproblem5(area,volume,jacflag,R0,H0); %find radius and height of a cylinder with given surface area and volume %note that there may be NO SOLUTION or MULTIPLE SOLUTIONS % %VARIABLES: %area: surface area of desired cylinder %volume: volume of desired cylinder %jacflag: set to 1 if want to calculate jacobian matrix explicitly, else 0 %R0,H0: initial guesses for radius and height %EXAMPLE: %consider a cylinder with R = 1 and H = 1... %area = 2*(pi*R^2) + (2*pi*R*H) = 4pi, volume = pi*R^2*H = pi; %SOUGHT SOLUTION %[R,H] = sampleproblem5(4*pi,pi,0,1.1,0.9) %[R,H] = sampleproblem5(4*pi,pi,1,1.1,0.9) %ANOTHER SOLUTION %[R,H] = sampleproblem5(4*pi,pi,1,0.5,3) %Initial guesses are important! %BEGIN CODE x0 = [R0,H0];%[volume^1/3,sqrt(area)]'; %initial guess tol = sqrt(eps)/100; %sqrt(eps) ~ 10^-8, often a reasonable tolerance param.volume = volume; param.area = area; format long; if(jacflag == 0) options = optimset('TolFun',tol,'TolX',tol,'Jacobian','off'); x = fsolve(@gcyl,x0,options,param); else options = optimset('TolFun',tol,'TolX',tol,'Jacobian','on'); x = fsolve(@gcyl_jac,x0,options,param); end R = x(1); H = x(2); return; function [y,Jacobian] = gcyl_jac(x,param); R = x(1); %cylinder radius H = x(2); %cylinder height y1 = 2*pi*R^2 + 2*pi*R*H - param.area; y2 = pi*R^2*H - param.volume; dy1dR = 4*pi*R + 2*pi*H; dy1dH = 2*pi*R; dy2dR = 2*pi*R*H; dy2dH = pi*R^2; Jacobian = [ dy1dR dy1dH ;... dy2dR dy2dH ]; y = [y1 y2]'; return; %if jacobian is set to off, then then the function is just: function y = gcyl(x,param); R = x(1); %cylinder radius H = x(2); %cylinder height y1 = 2*pi*R^2 + 2*pi*R*H - param.area; y2 = pi*R^2*H - param.volume; y = [y1 y2]'; return;