Determine Critical Points Numerically

I was wondering if anyone had any suggestions on how to find the critical points of this function numerically with respect to the variables r and phi. For the variable h, I do not have an exact value for it, but I know the range of values for it. I have been unable to solve this analytically using the solve function. Any help would be appreciated!

[Walter Roberson]

I tried copying the image into MATLAB to experiment with the function, but MATLAB doesn't seem to be able to execute images :(

[Ashmika Gupta]

EA=30000;
alpha=30;
beta=30;
gamma=180-(alpha+beta);
a=20;
n=6;
L_ab=a;
L_bc=(a*sind(alpha))/(sind(beta));
L_ac=(a*sind(alpha+beta))/(sind(beta));
hmax=L_bc*cosd(180-(90+(180-gamma)));

h = 0.1:0.1:hmax
syms r phi h
U=((n*EA)/2)*(L_ab*((2*(r/L_ab)*sin(pi/n)-1)^2)+L_bc*(((sind(beta)/sind(alpha))*sqrt(((h/L_ab)^2)...
    -(2*((r/L_ab)^2)*cos(phi))+2*((r/L_ab)^2)))-1)^2+L_ac*((sind(beta)/sind(alpha+beta))*sqrt(((h/L_ab)^2)...
    -(2*((r/L_ab)^2)*cos(phi+(2*pi/n)))+(2*((r/L_ab)^2)))-1)^2);

Would this help?

Answer 1

[Walter Roberson]

EA=30000;
alpha=30;
beta=30;
gamma=180-(alpha+beta);
a=20;
n=6;
L_ab=a;
L_bc=(a*sind(alpha))/(sind(beta));
L_ac=(a*sind(alpha+beta))/(sind(beta));
hmax=L_bc*cosd(180-(90+(180-gamma)));

%h = 0.1:0.1:hmax;

h = 0.1 + rand()*(hmax-0.1);
disp(h)

    1.3076

syms r phi
U=((n*EA)/2)*(L_ab*((2*(r/L_ab)*sin(pi/n)-1)^2)+L_bc*(((sind(beta)/sind(alpha))*sqrt(((h/L_ab)^2)...
    -(2*((r/L_ab)^2)*cos(phi))+2*((r/L_ab)^2)))-1)^2+L_ac*((sind(beta)/sind(alpha+beta))*sqrt(((h/L_ab)^2)...
    -(2*((r/L_ab)^2)*cos(phi+(2*pi/n)))+(2*((r/L_ab)^2)))-1)^2);

dUr = diff(U,r);

dUphi = diff(U,phi);

F = matlabFunction([dUr, dUphi], 'vars', {[r,phi]});
disp(F)

    @(in1)[in1(:,1).*9.0e+3+((sqrt(3.0).*sqrt(in1(:,1).^2.*cos(in1(:,2)+pi./3.0).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3))./3.0-1.0).*(in1(:,1)./1.0e+2-(in1(:,1).*cos(in1(:,2)+pi./3.0))./1.0e+2).*1.0./sqrt(in1(:,1).^2.*cos(in1(:,2)+pi./3.0).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3).*1.8e+6+(sqrt(in1(:,1).^2.*cos(in1(:,2)).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3)-1.0).*(in1(:,1)./1.0e+2-(in1(:,1).*cos(in1(:,2)))./1.0e+2).*1.0./sqrt(in1(:,1).^2.*cos(in1(:,2)).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3).*1.8e+6-1.8e+5,in1(:,1).^2.*sin(in1(:,2)+pi./3.0).*((sqrt(3.0).*sqrt(in1(:,1).^2.*cos(in1(:,2)+pi./3.0).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3))./3.0-1.0).*1.0./sqrt(in1(:,1).^2.*cos(in1(:,2)+pi./3.0).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3).*9.0e+3+in1(:,1).^2.*sin(in1(:,2)).*(sqrt(in1(:,1).^2.*cos(in1(:,2)).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3)-1.0).*1.0./sqrt(in1(:,1).^2.*cos(in1(:,2)).*(-1.0./2.0e+2)+in1(:,1).^2./2.0e+2+4.274324679548135e-3).*9.0e+3]

N = 100;
sols = zeros(N,4);

Nsol = 0;
opts = optimoptions(@fsolve, 'display', 'none');

for K = 1 : N
  guess = 10*randn(1,2);
  
  [onecrit, fval, exitflag] = fsolve(F, guess, opts);
  if exitflag > 0
      Nsol = Nsol + 1;
      sols(Nsol,:) = [guess, onecrit];
  end
end
sols = sols(1:Nsol,:);

disp(sols(1:10,:))

   -1.9004  -12.2623   -0.1945   -9.9571
   13.4320   -7.7723   20.0000  -11.5216
    8.7771  -15.8970   20.0000  -17.8048
   -0.3686    3.1735   -0.1945    2.6093
   -6.0614    2.3239   -8.4063    2.3963
    8.3212  -19.5582   29.1795  -19.2352
   -4.0384   -3.5463   -8.4063   -3.8869
    6.7465    9.1781   20.0000    7.3279
   -1.7631   12.8449   -0.1945   15.1756
   -2.9791   20.0151   -8.4063   21.2459

disp('---')

---

criticals = uniquetol(sols(:,3:4),0.0001,'byrows', true);
disp(size(criticals))

    28     2

format long g
disp(criticals)

         -8.40629280347899         -3.88687299162428
         -8.40629280347899         -10.1700582988039
         -8.40629280347899          2.39631231555531
         -8.40629280347899          8.67949762273489
         -8.40629280347899          21.2458682370941
        -0.194469915762064         -9.95709210493248
        -0.194469915762064         -3.67390679775289
        -0.194469915762064          2.60927850942669
        -0.194469915762064          15.1756491237859
          19.9780943191362         -7.33044737469427
          19.9780943191363         -1.04726206751469
          19.9780943191363           5.2359232396649
          19.9780943191363          -26.180003296233
          20.0000119522625          1.04472766456231
          20.0000119522625         -11.5216429497969
          20.0000119522625          7.32791297174189
          20.0000119522625          13.6110982789215
          20.0000119522625         -5.23845764261728
          20.0000119522625          19.8942835861011
          20.0000119522625         -17.8048282569765
          25.4510856003027         -12.5798701135576
          25.4510856003027          6.26968580798112
          25.4510856003027         -18.8630554207372
          25.4510856003027       -0.0134994991984623
           29.179493175519         -19.2351832373583
          29.1794931755192         -6.66881262299908
          29.1794931755192         -12.9519979301787
          29.1794931755192        -0.385627315819498