Solving Elliptic Integral Numerically
15 ビュー (過去 30 日間)
古いコメントを表示
Hello Everyone,
So I'm trying solve the following equations 10 and 15 from the paper called "Large deflection of a cantilever beam with multiple geometric and/or material discontinuities under a concentrated end‐point load". The equations can be seen below and I have uploaded a copy of the paper of as well.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/349034/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/349037/image.png)
I know what all the constants are (all the L's, E's, I's, F's, etc). I'm trying to find φ1 and φ2 and the way to go about this is to integrate the two equations (10 and 15), leave φ1 and φ2 in the resulting integrals, and since we know what L1 and L2 are, we can solve simultaneously. However due these functions being elliptic integrals and there is no closed form solution for these types of integrals. When I try to integrate the equation with L1 for example (not plugging in 0 and φ1), I receive the following (where p, p1, and p2 are φ, φ1, and φ2):
-(513^(1/2)*1000^(1/2)*((13*sin(p1) - 18*sin(p) + 5*sin(p2))/(13*sin(p1) + 5*sin(p2) - 18))^(1/2)*ellipticF(pi/4 - p/2, -36/(13*sin(p1) + 5*sin(p2) - 18)))/(500*(13*sin(p1) - 18*sin(p) + 5*sin(p2))^(1/2))
The "ellipticF" portion of the resulting integral does not allow me to solve for an actual result so I'm a little stuck with what I should do. I would appreciate any help or suggestions you could provide.
UPDATE
I have been working with the following script and for some reason I'm not getting the angles that are stated in the paper (φ1 and φ2 are .2935 and .5666). When I plug in the φ values from the paper into the L1 and L2 equations in script, I'm able to get the L1 and L2 values from the paper. However, when I simutaneously solve for both equations, I receive much larger values (please see in code). I would appreicate any help you could provide.
b=.0304;
F=4;
E1=200e9;
L1=.2;
h1=.001;
I1=b*h1^3/12;
E2=200e9;
L2=.1;
h2=.0005;
I2=b*h2^3/12;
A_L1=(2*E2*I2*F)/((E1*I1)^2);
B_L1=(2*F)/(E1*I1);
A_L2=sqrt(E2*I2/(2*F));
syms p1 p2 p
TSO_val = 5;
eq_L1 = sqrt(A_L1*(sin(p2)-sin(p1))+B_L1*(sin(p1)-sin(p)));
L1_int = int(1/eq_L1,p);
L1_int_p1 = subs(L1_int,p,p1);
L1_int_0 = subs(L1_int,p,0);
L1_int = L1_int_p1 - L1_int_0;
% L1_int1 = subs(L1_int, [p1, p2], [.2935, .5666]); % this confirms that my
% set up is right since when I solved this I get L1 = 20cm
% vpa(L1_int1)
eq_L2 = sqrt(sin(p2)-sin(p));
L2_int = A_L2*int(1/eq_L2,p);
L2_int_p2 = subs(L2_int,p,(p2-1e-30)); % I had to subtract by 1e-30 so didn't get a divison by 0 error
L2_int_p1 = subs(L2_int,p,p1);
L2_int = L2_int_p2 - L2_int_p1;
% L2_int1 = subs(L2_int, [p1, p2], [.2935, .5666]); % this confirms that my
% set up is right since when I solved this I get L2 = 10cm
% vpa(L2_int1)
[phi1,phi2] = vpasolve(L1_int == L1, L2_int == L2, [p1 p2])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% phi1 and phi2 end up being the following:
% phi1 =
%
% 51.385302992597424363840087698429 - 9.0368925680970086592118314144003i
%
%
% phi2 =
%
% 133.48826717725105017201958024896 - 9.2837426012587879832846341253788i
0 件のコメント
採用された回答
Alan Stevens
2020 年 8 月 22 日
The following gets close:
% Basic data
b=.0304;
F=4;
E1=200e9;
L1=.2;
h1=.001;
I1=b*h1^3/12;
E2=200e9;
L2=.1;
h2=.0005;
I2=b*h2^3/12;
% Constants used in integral calculations
K1 = sqrt(E2*I2/(2*F));
K2 = 2*E2*I2*F/(E1*I1)^2;
K3 = 2*F/(E1*I1);
K = [K1 K2 K3];
% Multi-dimensional Secant method
Phi = [0; pi]; % Initial guess
err = 1;
while err>10^-5
Phiold = Phi;
F = Ffn(Phi,K);
J = Jfn(Phi,K);
Phi = Phiold - J\F;
err = max(abs(Phi-Phiold));
end
disp(Phi)
% "Jacobian" function
function JJ = Jfn(Phi,K)
dphi = 10^-12;
Phi1 = Phi + dphi*[1;0];
Phi2 = Phi + dphi*[0;1];
JJ = [( Ffn(Phi1,K)-Ffn(Phi,K) )/dphi,...
( Ffn(Phi2,K)-Ffn(Phi,K) )/dphi];
end
% Function to be zeroed
function FF = Ffn(Phi,K)
L1=.2;
L2=.1;
phi1 = Phi(1);
phi2 = Phi(2);
K1 = K(1);
K2 = K(2);
K3 = K(3);
lo = 10^-16; % needed to protect against negative in square roots
factor1 = @(phi) max(K2*(sin(phi2)-sin(phi1)) + K3*(sin(phi1)-sin(phi)),lo);
factor2 = @(phi) max(sin(phi2)-sin(phi),lo);
kernel1 = @(phi) 1./sqrt( factor1(phi) );
kernel2 = @(phi) K1./sqrt( factor2(phi) );
FF = [integral(kernel1,0,phi1)-L1;
integral(kernel2,phi1,phi2)-L2];
end
producing Phi = [0.2936; 0.5669]
9 件のコメント
Walter Roberson
2023 年 8 月 21 日
Is the portion of the beam between s and L weightless? If not then there is load at L.
Prasanna Routray
2023 年 8 月 21 日
Ideally there is weight and it is not weightless between s and L. However, it confuses me to an extent that, I wonder how come the load is at one point and we assume it at the end.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460927/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1460932/image.jpeg)
The first figure shows that the load is at the end the maximum end-slope is at the end.
However, in the second figure, the load is at the midpoint and lets assume that as s<L. Here the maximum slopeoccurs at s=L/2 and remains same thereafter.
Cheers!
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Nonlinear Dynamics についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!