How to solve a system of 8 equations having integrals!

4 views (last 30 days)
Abubakar Rashid on 29 Nov 2022
Commented: Abubakar Rashid on 30 Nov 2022
d=1.18; E=27000000; L=121.53;
I=(pi*d^4)/64;
x1=16.99; x2=34.16; x3=51.33; x4=68.49; x5=85.66; x6=102.83; x7=115.61; x8=120.12;
F1=49; F2=49; F3=49; F4=49; F5=49; F6=49; F7=49; F8=11;
K1=100; K2=100; K3=100; K4=100; K5=100; K6=100; K7=100; K8=100;
These are my moment equation, that are further used in the deflection (d1, d2, ...., d8) calculations.
M1 = @(x)(F1*x-K1*d1*x);
M2 = @(x)(F1*x+F2*(x-x1)-K1*d1*x-K2*d2*(x-x1));
M3 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2));
M4 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3));
M5 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4));
M6 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)+F6*(x-x5)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4)-K6*d6*(x-x5));
M7 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)+F6*(x-x5)+F7*(x-x6)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4)-K6*d6*(x-x5)-K7*d7*(x-x6));
M8 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)+F6*(x-x5)+F7*(x-x6)+F8*(x-x7)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4)-K6*d6*(x-x5)-K7*d7*(x-x6)-K8*d8*(x-x7));
Below are my eight equations with defined boundaries for integrals. And it contains eight variables i.e. d1, d2, ...., d8 as well. I am fine with both numerical or analytical solution.
eq1 = d1 ==(1/(E*I))*(integral(@(x)(M1(x).*x),0,x1)+integral(@(x)(M2(x).*x),x1,x2)+integral(@(x)(M3(x).*x),x2,x3)+integral(@(x)(M4(x).*x),x3,x4)+integral(@(x)(M5(x).*x),x4,x5)+integral(@(x)(M6(x).*x),x5,x6)+integral(@(x)(M7(x).*x),x6,x7)+integral(@(x)(M8(x).*x),x7,x8));
eq2 = d2 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*(x-x1)),x1,x2)+integral(@(x)(M3(x).*(x-x1)),x2,x3)+integral(@(x)(M4(x).*(x-x1)),x3,x4)+integral(@(x)(M5(x).*(x-x1)),x4,x5)+integral(@(x)(M6(x).*(x-x1)),x5,x6)+integral(@(x)(M7(x).*(x-x1)),x6,x7)+integral(@(x)(M8(x).*(x-x1)),x7,x8));
eq3 = d3 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*(x-x2)),x2,x3)+integral(@(x)(M4(x).*(x-x2)),x3,x4)+integral(@(x)(M5(x).*(x-x2)),x4,x5)+integral(@(x)(M6(x).*(x-x2)),x5,x6)+integral(@(x)(M7(x).*(x-x2)),x6,x7)+integral(@(x)(M8(x).*(x-x2)),x7,x8));
eq4 = d4 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*(x-x3)),x3,x4)+integral(@(x)(M5(x).*(x-x3)),x4,x5)+integral(@(x)(M6(x).*(x-x3)),x5,x6)+integral(@(x)(M7(x).*(x-x3)),x6,x7)+integral(@(x)(M8(x).*(x-x2)),x7,x8));
eq5 = d5 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*(x-x4)),x4,x5)+integral(@(x)(M6(x).*(x-x4)),x5,x6)+integral(@(x)(M7(x).*(x-x4)),x6,x7)+integral(@(x)(M8(x).*(x-x4)),x7,x8));
eq6 = d6 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*0),x4,x5)+integral(@(x)(M6(x).*(x-x5)),x5,x6)+integral(@(x)(M7(x).*(x-x5)),x6,x7)+integral(@(x)(M8(x).*(x-x5)),x7,x8));
eq7 = d7 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*0),x4,x5)+integral(@(x)(M6(x).*0),x5,x6)+integral(@(x)(M7(x).*(x-x6)),x6,x7)+integral(@(x)(M8(x).*(x-x6)),x7,x8));
eq8 = d8 ==(1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*0),x4,x5)+integral(@(x)(M6(x).*0),x5,x6)+integral(@(x)(M7(x).*0),x6,x7)+integral(@(x)(M8(x).*(x-x7)),x7,x8));

Torsten on 29 Nov 2022
d0 = ones(1,8);
d = fsolve(@fun,d0)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
d = 1×8
0.5346 0.5133 0.4823 0.4784 0.3431 0.2181 0.0778 0.0067
norm(fun(d))
ans = 1.4844e-15
function res = fun(D)
d=1.18;
E=27000000;
L=121.53;
I=(pi*d^4)/64;
x1=16.99; x2=34.16; x3=51.33; x4=68.49; x5=85.66; x6=102.83; x7=115.61; x8=120.12;
F1=49; F2=49; F3=49; F4=49; F5=49; F6=49; F7=49; F8=11;
K1=100; K2=100; K3=100; K4=100; K5=100; K6=100; K7=100; K8=100;
d1 = D(1);
d2 = D(2);
d3 = D(3);
d4 = D(4);
d5 = D(5);
d6 = D(6);
d7 = D(7);
d8 = D(8);
M1 = @(x)(F1*x-K1*d1*x);
M2 = @(x)(F1*x+F2*(x-x1)-K1*d1*x-K2*d2*(x-x1));
M3 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2));
M4 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3));
M5 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4));
M6 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)+F6*(x-x5)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4)-K6*d6*(x-x5));
M7 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)+F6*(x-x5)+F7*(x-x6)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4)-K6*d6*(x-x5)-K7*d7*(x-x6));
M8 = @(x)(F1*x+F2*(x-x1)+F3*(x-x2)+F4*(x-x3)+F5*(x-x4)+F6*(x-x5)+F7*(x-x6)+F8*(x-x7)-K1*d1*x-K2*d2*(x-x1)-K3*d3*(x-x2)-K4*d4*(x-x3)-K5*d5*(x-x4)-K6*d6*(x-x5)-K7*d7*(x-x6)-K8*d8*(x-x7));
res(1) = d1 -((1/(E*I))*(integral(@(x)(M1(x).*x),0,x1)+integral(@(x)(M2(x).*x),x1,x2)+integral(@(x)(M3(x).*x),x2,x3)+integral(@(x)(M4(x).*x),x3,x4)+integral(@(x)(M5(x).*x),x4,x5)+integral(@(x)(M6(x).*x),x5,x6)+integral(@(x)(M7(x).*x),x6,x7)+integral(@(x)(M8(x).*x),x7,x8)));
res(2) = d2 -((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*(x-x1)),x1,x2)+integral(@(x)(M3(x).*(x-x1)),x2,x3)+integral(@(x)(M4(x).*(x-x1)),x3,x4)+integral(@(x)(M5(x).*(x-x1)),x4,x5)+integral(@(x)(M6(x).*(x-x1)),x5,x6)+integral(@(x)(M7(x).*(x-x1)),x6,x7)+integral(@(x)(M8(x).*(x-x1)),x7,x8)));
res(3) = d3 -((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*(x-x2)),x2,x3)+integral(@(x)(M4(x).*(x-x2)),x3,x4)+integral(@(x)(M5(x).*(x-x2)),x4,x5)+integral(@(x)(M6(x).*(x-x2)),x5,x6)+integral(@(x)(M7(x).*(x-x2)),x6,x7)+integral(@(x)(M8(x).*(x-x2)),x7,x8)));
res(4) = d4 -((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*(x-x3)),x3,x4)+integral(@(x)(M5(x).*(x-x3)),x4,x5)+integral(@(x)(M6(x).*(x-x3)),x5,x6)+integral(@(x)(M7(x).*(x-x3)),x6,x7)+integral(@(x)(M8(x).*(x-x2)),x7,x8)));
res(5) = d5-((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*(x-x4)),x4,x5)+integral(@(x)(M6(x).*(x-x4)),x5,x6)+integral(@(x)(M7(x).*(x-x4)),x6,x7)+integral(@(x)(M8(x).*(x-x4)),x7,x8)));
res(6) = d6 -((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*0),x4,x5)+integral(@(x)(M6(x).*(x-x5)),x5,x6)+integral(@(x)(M7(x).*(x-x5)),x6,x7)+integral(@(x)(M8(x).*(x-x5)),x7,x8)));
res(7) = d7 -((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*0),x4,x5)+integral(@(x)(M6(x).*0),x5,x6)+integral(@(x)(M7(x).*(x-x6)),x6,x7)+integral(@(x)(M8(x).*(x-x6)),x7,x8)));
res(8) = d8 -((1/(E*I))*(integral(@(x)(M1(x).*0),0,x1)+integral(@(x)(M2(x).*0),x1,x2)+integral(@(x)(M3(x).*0),x2,x3)+integral(@(x)(M4(x).*0),x3,x4)+integral(@(x)(M5(x).*0),x4,x5)+integral(@(x)(M6(x).*0),x5,x6)+integral(@(x)(M7(x).*0),x6,x7)+integral(@(x)(M8(x).*(x-x7)),x7,x8)));
end
Abubakar Rashid on 30 Nov 2022
Thank you Torsten, it helps a lot!

Walter Roberson on 29 Nov 2022
Do not use == with integral() . integral() is for numeric integration, and if you use numeric integration returning a numeric value then the right side of the == would be numeric, and the == would be for asking whether d1 is exactly the same value, bit-for-bit identical (except for different nans being considered the same, and negative 0 being treated the same as regular 0)
If you want to process this system numerically using fsolve() then you should change the NAME == EXPRESSION into NAME - (EXPRESSION) and then you would have the function return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8]
You could alternately switch to using the Symbolic Toolbox and use int() or vpaintegral() instead of integral, and try to solve() or vpasolve() the vector of equations (which can be written with == in that case.)
Abubakar Rashid on 30 Nov 2022
Thank you Walter for the quick feedback.

Categories

Find more on Calculus in Help Center and File Exchange

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by