Solving Elliptic Integral Numerically

Question

Uche Agwu 2020-8-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/582191-solving-elliptic-integral-numerically

评论： Prasanna Routray 2023-8-21

采纳的回答： Alan Stevens

Ramezanpour-Nahvi2017_Article_LargeDeflectionOfACantileverBe.pdf

在 MATLAB Online 中打开

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.

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2020-8-22

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/582191-solving-elliptic-integral-numerically#answer_483461

在 MATLAB Online 中打开

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 个评论
显示 7更早的评论隐藏 7更早的评论

Alan Stevens 2020-8-22

移动：Walter Roberson 2023-8-21

在 MATLAB Online 中打开

The method is basically a slight modification of the Newton-Raphson method for finding algebraic roots. In one dimension Newton-Raphson iterates towards a solution of f(x) = 0, using

x(n+1) = x(n) - f(x(n))/dfdx(x(n))

This can only be done when dfdx(x) is known. If there is no analytical (or easily obtainable) expression for dfdx, then a discrete approximation can be used, namely:

dfdx = (f(x+deltax) - f(x))/deltax

This is then known as the secant method (because one is approximating the gradient of a curve by a secant).

This is what I did for your case, except I used a multi-dimensional version, in which dfdx is given by the "Jacobian" (strictly, the Jacobian would involve the true derivatives). You had two functions say, f1 and f2, that had to be zeroed by finding the appropriate values of x1 and x2. The true Jacobian would be a matrix J = [df1/dx1 df1/dx2; df2/dx1 df2/dx2], and the iteration is x(n+1) = x(n) - J^-1*f, where x and f are 2x1 vectors, and J is a 2x2 matrix. In our case, the dfdx's are approximated by expressions like (f(x+deltax) - f(x))/deltax (for x read phi, for deltax read dphi) (in Matlab J^-1*f is best calculated using the backslash operator: J\f).

I didn't spend any time trying to debug your program as I couldn't see the point of using symbolic operations here to obtain a purely numeric solution.

Prasanna Routray 2023-8-10

移动：Walter Roberson 2023-8-21

在 MATLAB Online 中打开

Hi Alan,

I ran your code and it ran fine.

I did the modifications according to the elliptical equation (10) given in https://iopscience.iop.org/article/10.1088/0143-0807/23/3/317

Trying to solve for phi0 for a given value of E, I, F, and L

The adapted code goes like this:

% Basic data
b=.0304;
F=4;
E1=200e9;
L1=.5;
h1=.001;
I1=b*h1^3/12;
% Constants used in integral calculations
K1 = sqrt(E2*I2/(2*F));
% Multi-dimensional Secant method 
Phi = 0; % Initial guess
err = 1;
while err>10^-5
    Phiold = Phi;
    F = Ffn(Phi,K1);
    J = Jfn(Phi,K1);
    Phi = Phiold - J\F;
    err = max(abs(Phi-Phiold));
end
disp(Phi)
% "Jacobian" function
function JJ = Jfn(Phi,K1)
dphi = 10^-12;
Phi0 = Phi + dphi;
JJ = ( Ffn(Phi0,K1)-Ffn(Phi,K1) )/dphi;
end
% Function to be zeroed
function FF = Ffn(Phi,K1)
L1=.2;
phi0 = Phi(1);
lo = 10^-16;  % needed to protect against negative in square roots
factor1 = @(phi) max((sin(phi0)-sin(phi)),lo);
kernel1 = @(phi) K1./sqrt( factor1(phi) );
FF = integral(kernel1,0,phi0)-L1;
end

I did not understand much but just tried to run it for a special case for which you have already given a code.

Could you please give some feedback so that I can see it working.

Cheers!

Prasanna

Walter Roberson 2023-8-21

@Prasanna Routray

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.

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!

请先登录，再进行评论。

Solving Elliptic Integral Numerically

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving Elliptic Integral Numerically

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论