Solving Elliptic Integral Numerically

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

Answer 1

[Alan Stevens]

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]

[Uche Agwu]

Oh interesting, thanks this is definitely an interesting way to do it! Can you briefly explain your method/thought process? How did you know to use the secant method and creating a jacobian function? Also what is "dphi" and why do you divide the functions in "JJ" by "dphi"? Finally do you know why my script wasn't working or maybe not the best method?

[Alan Stevens]

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.

[Uche Agwu]

Oh I see that makes sense now. Thanks so much for the explanation and I truly appreciate all the help!

[Prasanna Routray]

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

[Alan Stevens]

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(E1*I1/(2*F));   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multi-dimensional Secant method 
Phi = pi/3; % 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)

    0.1565

% Derivative 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

[Prasanna Routray]

Hi Alan, I could solve that one using a different initial guess. I got the same result. Thanks for your valuable contribution and I really appreciate the effort you put for quick and clear response.

Cheers!

[Prasanna Routray]

Hi Alan,

I'm back here with a silly question. This is not related to matlab but related to the mecanics of beam itself.

I should not be asking MATLAB irrelevant questions but I believe this would be a quicker way to get a solution to my doubt.

The literatures shared in this thread assume that we have load at the endpoint (so the effective length is L). How about the case when the load is not at the end-point and is anywhere along the beam (s<=L).

How do I modify the expression to calculate Phi0?

What I did is I used equation (9) instaed of equation (10) and assume that the load is applied at 's' and not 'L'.

I do not have a clear intution about the feasibility of this to calculate Phi0.

Cheers!

[Walter Roberson]

@Prasanna Routray

Is the portion of the beam between s and L weightless? If not then there is load at L.

[Prasanna Routray]

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!