Solving a system of nonlinear equations

Im working on solving for the static equilibrium of 3 connected hanging masses. This is shown in the diagram.

The masses and lengths of string are known. I end up with a system of nonlinear equations. I have 8 unknowns to solve for, which are 4 angles (t1, t2, t3, t4) and 4 tensions (T1, T2, T3, T4).

So far i have written the following code:

clear

m=[1 1 1]; %arbitrary values for mass

l=[5 5 5 5]; %arbitrary values for rope length

g=9.81;

D=20; %length of bar

syms t1 t2 t3 t4 T1 T2 T3 T4

x= [ sin(t1) sin(t2) sin(t3) sin(t4) cos(t1) cos(t2) cos(t3) cos(t4) T1 T2 T3 T4]';

%the following are the 12 equations that can be obtained from the system

f1= l(1).*x(5) + l(2).*x(6)+l(3).*x(7)+l(4).*x(8)-20;

f2= l(1).*x(1) + l(2).*x(2) - l(3).*x(3) - l(4).*x(4);

f3= -x(9).*x(5) + x(10).*x(6);

f4= m(1).*g + x(10).*x(2) - x(9).*x(1);

f5= x(11).*x(7) - x(10).*x(6);

f6= m(2).*g -x(10).*x(2)-x(11).*x(3);

f7= x(12).*x(8) -x(11).*x(7);

f8= m(3).*g + x(11).*x(3) - x(12).*x(4);

f9= (x(1)).^(2) + (x(5)).^(2) -1;

f10= (x(2)).^(2) + (x(6)).^(2) -1;

f11= (x(3)).^(2) + (x(7)).^(2) -1;

f12= (x(4)).^(2) + (x(8)).^(2) -1;

f = [f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12];

-------------------------------------------------

now i understand that fsolve cannot be used for symbolic variables, so how do i go about solving f?

[Zayan Ahsan Ali]

Update: I've been able to solve it (i think) with vpasolve:

[sol_t1, sol_t2, sol_t3, sol_t4, sol_T1, sol_T2, sol_T3, sol_T4] = vpasolve(f, [t1 t2 t3 t4 T1 T2 T3 T4])

The solutions i get for tensions are fine, but the angles are strange, some of which are very large. I'm guessing that theres a pi-integer difference, and i'm getting solutions in the wrong region. How do i change that to get solutions within 0 to 90 degree?

[Walter Roberson]

You can add a third parameter to vpasolve() which is a numeric array with two columns, with one row for each variable (in the same order being solved for.) The first column is lower bounds; the second column is upper bounds. Use nan in both columns for entries that should not be bounded.

Caution: vpasolve() uses a Newton-raphson type solver, which overshoots by necessity. If an overshoot on a solution near (but within) a boundary goes out of bounds, vpasolve() will stop solving there rather than moving the value to inside the boundary.

[Walter Roberson]

Note:

You have x = [....]' which is conjugate transpose. You do not make use of the fact that x is a column vector, and you appear to be wanting to work with real values, suggesting that you meant non-conjugate transpose at most. If you remove the conjugate part or declare the symbols to be real valued, then the last four entries of f vanish to 0. Your f is 12 long, of which 8 are non-zero, so you should probably get rid of those last 4.

[Walter Roberson]

Odd, testing in a different software package, I find it claiming that there are a series of 480 closed form complex-valued solutions based upon roots of polynomials, the main one of which is degree 20. And every single value involved in those 480 solutions is complex valued: not even one real-valued solutions.

The numeric values you get with normal vpasolve() for T1, T2, T3, T4 are sort of large magnitude; if you subs() the solutions back into f, you will see that some of the results are not zero, partly due to round off error acting combining with large numeric values. But the problems persist even if you move to digits(100): f(4) and f(8) are particular problems.

[Walter Roberson]

As I go for higher precision such as 250 digits, the error is not vanishing for f(4) and f(8), and can even be the opposite sign from before.

It should probably be assumed that vpasolve() is returning false roots.

[John D'Errico]

Hmm. I thought clouds were really light? They float, don't they? ;-)

Seriously, it seems interesting in that these problems are usually formulated as a truss, with pinned joints, but bars that can vary in length. So just a force balance computation, basic finite elements. What you might see in a first course in FEM. It allows you to solve the problem using a purely linear system of equations, in terms of the coordinates of the three pinned joints.

Here, you are assuming bars of fixed lengths, so essentially infinite stiffness, an infinite Young's modulus. I have a funny feeling that you could formulate the problem in the form of a classical truss where the bars have a finite Young's modulus, and then just obtain the desired solution as the limit as E-->inf. The two methodologies should yield the same result in the end, but the scheme I described will not require you to solve a nonlinear system of equations.

Of course at the end, once you have the final set of displacements, you can back out the tensile stresses in each bar, so only the displacements for each node are necessary to describe the state of the system.

[Alex Sha]

An approximate solution：

x1: -0.000386859711469918
x2: -0.000128952090745106
x3: -0.000128955577238622
x4: -0.000386863220600714
x5: 0.999999947230989
x6: 1.00000012310153
x7: 0.999999962070972
x8: 0.999999952070866
x9: -38036.8719191646
x10: -38036.865229583
x11: -38036.8713547058
x12: -38036.8717350701

Fevl:
-7.76282149672625E-8
3.49781215599771E-8
1.71785359270871E-8
8.78310757457257E-10
-2.49274307861924E-8
1.35815358959235E-9
8.46193870529532E-9
8.91805740366181E-10
4.41224170533161E-8
2.6283171705721E-7
-5.92285137601678E-8
5.38048858800266E-8

Answer 1

[John D'Errico]

he he. I just noticed a major flaw with what you are trying to compute.

l=[5 5 5 5]; %arbitrary values for rope length
g=9.81;
D=20; %length of bar

You are formulating the problem as a set of 4 bars, with fixed lengths, l1=l2=l3=l4=5.

They hang from a pair of pinned joints, set 20 units apart. Now, you are asking to compute the angles, based on bars/rope with essentially infinite stiffness, so no extension is possible. Note that sum(l)==D.

The only real solution that exists to the problem you are trying to solve is the one with all 4 members in an exact line, thus theta1=theta2=theta3=theta4=0. (The stresses in each bar/rope will be interesting, of course. I wonder what they are made of, since it will be a strong enough material you could make a space elevator from it. None of that weak carbon nanotube crap here. My guess is those members might be made from unobtainium.)

The point is, if you do try to solve the problem as you formulated it, you will probably get garbage results, thus perhaps why you are having problems. Instead, you may need to use some reasonable bar lengths instead, perhaps [6 6 6 6] might be a better idea.

Always remember that a solver will try to compute what you tell it to compute. It won't worry that you are trying to compute something physically meaningless. You are the one who needs to think about the physical meaning of the results, to verify they actually correspond to something realistic.

Even with all that, I might wonder if the better approach is as I suggested in my comment. Thus, formulate this as a simple force balance, with an unknown Young's modulus E for the rope stiffnesses. Then look at the solution as E-->inf. The system should be a linear system of equations, solvable for an unknown modulus. Then look at the limit.

[Zayan Ahsan Ali]

Yes, i noticed that error aswell, a bit silly really. Better input values seem to have fixed my tension values :) Thanks alot!

[John D'Errico]

The most difficult problems can sometimes have simple resolutions. :)