Determine Critical Points Numerically

2 次查看（过去 30 天）
Ashmika Gupta
Ashmika Gupta 2021-2-23
回答： Walter Roberson 2021-2-23
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!
  2 个评论
Walter Roberson
Walter Roberson 2021-2-23
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
Ashmika Gupta 2021-2-23
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2021-2-23
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by