Computing the explicit expression of the integral of a function of two variables inetgrated with respective to only one variable

Hello Madam or Sir,

I am computing the Lagrange equation that models the dynamic behavior of flexible link manipulator, in my code below I am having a issue for getting the explicit expression of the Kinetic energy. The inetgral of the function of two variables but inetgrated with respective to only one variable does not give me the explicit answer.

Here is my code:

%% 1. Parameters 
syms L1 m1 J1 rho E I k_stif1 L0 theta1 theta2 real
syms L2 m2 J2 rho E k_stif2 real
syms mp JhP real
syms Jh1 mh1 Jh2 mh2 real

syms eta_m eta_g km kg Rm b1 b2 Vm1 Vm2 real
syms q1 q2 q3 q4 q5 q6 q1d q2d q3d q4d q5d q6d q1dd q2dd q3dd q4dd q5dd q6dd x t real
syms x1(t) x1d x1dd x2(t) x2d x2dd x3d x3dd x4d x4dd;
syms Q11(t) Q12(t) Q21(t) Q22(t);
x1d  = diff(x1(t),t);     x2d = diff(x2(t),t);
x1dd = diff(x1(t),t,t);  x2dd = diff(x2(t),t,t);
x3d  = diff(Q11(t),t);     x4d = diff(Q12(t),t);
x3dd = diff(Q11(t),t,t);  x4dd = diff(Q12(t),t,t);
x5d  = diff(Q21(t),t);     x6d = diff(Q22(t),t);
x5dd = diff(Q21(t),t,t);  x6dd = diff(Q22(t),t,t);

%% 2. Numerical values
m1 = 0.065; m2 = 0.070; mh1 = 0.08; mh2 = 0.05; L0 = 0.0;   L1 = 0.3;   L2 = 0.2;   J1 = 0.00195;                                                            
J2 = 0.000933;  Jh1 = 0.099; Jh2 = 0.092;  fc = 3.2;  k_stif1 = 2*pi*fc;  k_stif2 = 2*pi*fc;    
eta_m = 0.69;   eta_g = 0.9;    km = 0.00767;   kg = 70;    Rm = 2.6;   b1 = 1.99;  b2 = 1.99;
mP = 0.69; JhP = 2.08e-3;   E = 1;  I = 1;  rho = 0.2;  be1 = [6.24 15.6]; be2 = [15.2 42.0];

%% 2. Vectorize symbolic variables
L = [L0;L1;L2];                                                                                           % vector length of links
J = [J1;J2];                                
k = [k_stif1;k_stif2];
m = [m1;m2;mP];
Jh = [Jh1;Jh2;JhP];
mh = [mh1;mh2];
q = [q1;q2;q3;q4;q5;q6];
qdot = [q1d;q2d;q3d;q4d;q5d;q6d];
qddot = [q1dd;q2dd;q3dd;q4dd;q5dd;q6dd];
X = [x1(t);x2(t);Q11(t);Q12(t);Q21(t);Q22(t)];
Xdot = diff(X,t);
Xddot = diff(Xdot,t);
u = [Vm1;Vm2];
b = [b1;b2];
omega1 = zeros(length(be1),1);
omega2 = zeros(length(be2),1);
g1 = zeros(length(be1),1);
g2 = zeros(length(be2),1);
Ybar1=sym(zeros(length(be1),1)); C1=zeros(length(be1),1);
Ybar2=sym(zeros(length(be2),1)); C2=zeros(length(be2),1);
3.Input parameters
Tau1 = ((eta_m*eta_g*km*kg/Rm)*(u(1)-km*kg*diff(X(1),t)))-b(1).*diff(X(1),t);
Tau1 = subs(Tau1,Xdot,qdot);
Tau2 = ((eta_m*eta_g*km*kg/Rm)*(u(2)-km*kg*diff(X(2),t)))-b(2).*diff(X(2),t);
Tau2 = subs(Tau2,Xdot,qdot);
Tau3 = 0;
Tau4 = 0;
Tau5 = 0;
Tau6 = 0;
f = vpa([Tau1;Tau2;Tau3;Tau4;Tau5;Tau6],3)

%% 4. Deflection y(x,t)
for j=1:length(be1)
    omega1(j) = (be1(j)/L(2)).^2.*sqrt(E*I/rho);
    g1(j) = vpa(((sin(be1(j)) - sinh(be1(j))) + (m(1)*be1(1)/(rho*L(2)))*(cos(be1(j)) - cosh(be1(j))))/((cos(be1(j)) + cosh(be1(j))) - (m(1)*be1(j)/(rho*L(2))).*(sin(be1(j)) - sinh(be1(j)))),3);
    Ybar1(j) = vpa(cos(be1(j)*x/L(2)) - cosh(be1(j)*x/L(2)) + g1(j)*(sin(be1(j)*x/L(2)) - sinh(be1(j)*x/L(2))));
    C1(j) = vpa(1/sqrt(int(rho*Ybar1(j)^2,x,[0 L(2)]) + m(1)*compose(Ybar1(j)^2,L(2))));
end
y1 = vpa(C1(1).*Ybar1(1).*X(3) + C1(2).*Ybar1(2).*X(4))


for j=1:length(be2)
    omega2(j) = (be2(j)/L(3)).^2.*sqrt(E*I/rho);
    g2(j) = vpa(((sin(be2(j)) - sinh(be2(j))) + (m(2)*be2(1)/(rho*L(3)))*(cos(be2(j)) - cosh(be2(j))))/((cos(be2(j)) + cosh(be2(j))) - (m(2)*be2(j)/(rho*L(3))).*(sin(be2(j)) - sinh(be2(j)))),3);
    Ybar2(j) = vpa(cos(be2(j)*x/L(3)) - cosh(be2(j)*x/L(3)) + g2(j)*(sin(be2(j)*x/L(3)) - sinh(be2(j)*x/L(3))));
    C2(j) = vpa(1/sqrt(int(rho*Ybar2(j)^2,x,[0 L(3)]) + m(2)*compose(Ybar2(j)^2,L(3))));
end
y2 = vpa(C2(1).*Ybar2(1).*X(5) + C2(2).*Ybar2(2).*X(6));

%% 5. Transformation matrices
A1 = [cos(x1(t)) -sin(x1(t)); sin(x1(t)) cos(x1(t))];               %Rigide Rotation Matrix Transformation from (X0,Y0) to (X1,Y1)
A2 = [cos(x2(t)) -sin(x2(t)); sin(x2(t)) cos(x2(t))];               %Rigide Rotation Matrix Transformation from (X0,Y0) to (X1,Y1)
E0 = [1, 0; 0, 1];                                                  %rotation matrix of the flexible link 1 at the end point
E1 = [1, -subs(diff(y1,x),x,L(2)); subs(diff(y1,x),x,L(2)), 1];      %rotation matrix of the flexible link 2 at the end point
T0 = eye(2);
T1 = T0*E0*A1;
T2 = T1*E1*A2;

%% 6. Links and Tipload positions
p01 = subs([x ; y1],x,L(1));
r11 = [x ; y1];
r01 = p01 + T1*r11;
p12 = subs([x ; y1],x,L(2));
p02 = p01 + T1*p12; 
r22 = [x ; y2];
r02 = p02 + T2*r22;
p23 = subs([x ; y2],x,L(3));
p03 = p02 + T2*p23;
r03 = p03;
alpha1 = X(1) + subs(diff(y1,x),x,L(1)); 
alpha2 = X(1) + X(2) + subs(diff(y1,x),x,L(2));
alpha3 = X(1) + X(2) + subs(diff(y1,x),x,L(2)) + subs(diff(y2,x),x,L(3))

%% 7. Angular velocity of links and tipload
W1 = diff(alpha1,t);
W2 = diff(alpha2,t)
W3 = diff(alpha3,t)
V1 = diff(r01,t)
V2 = diff(r02,t)
VM = diff(r03,t)

%% 8. Kinetic Energy of Hub, Links and Tipload
TH1 = (1/2)*Jh(1)*W1^2 + (1/2) * ((diff(p01,t))') * mh(1) * (diff(p01,t));
TH2 = (1/2)*Jh(2)*W2^2 + (1/2) * ((diff(p02,t))') * mh(2) * (diff(p02,t))
TL1 = (1/2)*int(rho * (V1') * V1,x,0,L(2))
TL2 = (1/2)*int(rho * (V2') * V2,x,0,L(3)) 
TP = (1/2)*Jh(3)*W3^2 + (1/2) * ((diff(p03,t))') * m(3) * (diff(p03,t))
K = TH1 + TH2 + TL1 + TL2 + TP;

When computing the code, from the result of TL1 is shown as :

I am wonderng why Matlab can not compute the explicit answer.

Any help please ?

Answer 1

[Walter Roberson]

f = vpa([Tau1;Tau2;Tau3;Tau4;Tau5;Tau6],3)

You are using finite precision software floating point numbers.

TL1 = (1/2)*int(rho * (V1') * V1,x,0,L(2))

You are asking for indefinitely precise closed-form solutions for expressions that are constructed from finite precision software floating point numbers.

m1 = 0.065;

Written as a formula that represents the set of numbers that round to 65/1000 -- the set of numbers 65/1000 - 5/10000 inclusive to 65/1000 + 5/10000 exclusive. You have two choices when you work symbolically with such numbers:

1. Convert to a rational number with a uncertainty variable, such as syms delta_m1; assume(-5/10000 <= delta_m1 & delta_m1 < 5/10000); m1 = 65/1000 + delta_m1; Or
2. Do not use int() or solve() or symsum() on any expression involving those numbers, as those functions are logically only for indefinite precision closed form solutions

TL1 = (1/2)*int(rho * (V1') * V1,x,0,L(2))

Do you definitely intend complex conjugate transpose for V1, or do you just intend transpose ?

[Cianyi Kazadi]

TL1 = (1/2)*int(rho * (V1') * V1,x,0,L(2))

here I just intend transpose opereator.

f = vpa([Tau1;Tau2;Tau3;Tau4;Tau5;Tau6],3)

At this stage, f does not affect the computation of TL1 since f appears only at the end of the Lagrange equation.

2. Do not use int() or solve() or symsum() on any expression involving those numbers, as those functions are logically only for indefinite precision closed form solutions

I do not get this point, you intend that int() or solve() can not give me an exlicit expression ?

Thanks for you reply

[Walter Roberson]

(rho * (V1')) * V1 involves the undefined functions Q11(t), Q12(t), x1(t)

(rho * (V2')) * V2 involves the undefined functions Q11(t), Q12(t), Q21(t), Q22(t), x1(t), x2(t)

You cannot expect a closed form solution for an integral of an expression that involves an unknown function.

You introduce those functions into V1 and V2 by the way you define y1 and y2 by multiplying terms by various X() values, where X has been defined as a vector of derivatives.

I see no reason to hope that the situation could ever be salvaged.