How to solve a system of 8 equations having integrals!

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));

Answer 1

[Torsten]

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]

Thank you Torsten, I appreciate your effort.

Can you please explain briefly for my better understanding. Thank you!

And I am assuming this is an numerical solution using fsolve. Right?

[Torsten]

The code solves your equations eq1,...,eq8 for d1,...,d8 using Newton's method (thus a numerical solution).

I think equations 1-8 are just linear equations in the unknowns d1,...,d8. Thus in principle you solve a linear system of equations of the form A*d = b.

But getting the entries of A and b from your equations would be quite tedious - so I solved it as if it were a nonlinear system.

A symbolic approach and EquationsToMatrix could be of help to recover A and b.

[Abubakar Rashid]

Thank you Torsten, it helps a lot!

Answer 2

[Walter Roberson]

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.)

[Walter Roberson]

I recommend you get rid of those terms such as integral(@(x)(M1(x).*0),0,x1) . You know that the answer has to be scalar 0 -- if x1 were a vector then integral() would complain about the bounds needing to be a scalar, and if M1 returned a vector for a scalar input then integral() would complain about the function returning the wrong size. Don't make the function do more work than it needs to do.

[Abubakar Rashid]

Thank you Walter for the response.

Can you suggest me a single approach to follow, as I am getting confused. Pardon me for my lil experience with Matlab.

[Walter Roberson]

Torsten's response is good, but for efficiency I would change, for example,

  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)));

to

   res(2) = d2 -((1/(E*I))*(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)));
 

which got rid of the integral(@(x)(M1(x).*0),0,x1) since we know the result is going to be 0 for the term.

[Abubakar Rashid]

Thank you Walter for the quick feedback.