# How do I solve vector loop equations

This is a simple case, representing a slider mechanism. Crank length Ra, connecting rod length Rb. Crank angle is given as ta. Then gudgen pin displacement from crank origin = x and conrod angle = tb

Resolving in x and one in y dirextions we have two equations in two unknowns. These should be simple to solve!

I tried the following on sagecell.sagemath.org:

Ra=0.5; Rb=3.0; ta=pi/2; var('x tb');
sol = solve([Ra*cos(ta)+Rb*cos(tb)==x, Ra*sin(ta)+Rb*sin(tb)==0.0],x,tb,solution_dict=True)
sol

[{3.0*cos(tb): x}, {3.0*sin(tb) + 0.5: 0.0}]


The result is obviously nonsense.

What am I doing wrong?

edit retag close merge delete

Sorry I don't seem able to format the paste from save

( 2023-01-13 21:39:47 +0100 )edit

The code is the following (brad.can, please confirm)

Ra=0.5;
Rb=3.0;
ta=pi/2;
var('x tb');
sol = solve([Ra*cos(ta)+Rb*cos(tb)==x, Ra*sin(ta)+Rb*sin(tb)==0.0],x,tb,solution_dict=True)
sol

( 2023-01-14 06:49:36 +0100 )edit

Yes correct. But the sage result is:

[{3.0cos(tb): x}, {3.0sin(tb) + 0.5: 0.0}]

The numerical result should be x=2.958 and tb = -0.167 (MathCAD solutions) If I substitute tb 3.0cos(-0.167) is indeed 2.958 and 3.0sin(-0.167)+0.5 is 0.0

So we seem to have a partial result in the dict but key and value of the first is swopped and the second says zero equals zero!

It's clear that solve does not like these particular equations perhaps because the second could just be rearanged for tb. So in this case we don't need solve at all.

( 2023-01-14 11:24:26 +0100 )edit

The result is oviously nonsense.

No. the result is a re-statement of your original system, meaning that Sage's default solver (i.e. Maxima's) can't solve it.

But there are solutions : see my answer below...

( 2023-01-14 11:51:04 +0100 )edit

Sort by ยป oldest newest most voted

Well...

sage: var('x tb');
sage: Ra=1/2; Rb=3; ta=pi/2;
sage: Sys = [Ra*cos(ta)+Rb*cos(tb)==x, Ra*sin(ta)+Rb*sin(tb)==0] ; Sys
[3*cos(tb) == x, 3*sin(tb) + 1/2 == 0]


Sage's default solver (i. e. Maxima's) isn't especially good at solving trigonometric systems :

sage: solve(Sys, [x, tb])
[3*cos(tb) == x, 3*sin(tb) + 1/2 == 0]
sage: solve(Sys, [x, tb], to_poly_solve=True)
[3*cos(tb) == x, 3*sin(tb) + 1/2 == 0]
sage: solve(Sys, [x, tb], to_poly_solve="force")
[3*cos(tb) == x, 3*sin(tb) + 1/2 == 0]


but Sympy's can solve this one :

sage: solve(Sys, [x, tb], algorithm="sympy")
[{tb: pi + arcsin(1/6), x: -1/2*sqrt(35)}, {tb: -arcsin(1/6), x: 1/2*sqrt(35)}]


(and, FWIW, so can the gratis Wolfram engine :

sage: mathematica.Solve(Sys, [x, tb])
{{x -> -Sqrt[35]/2, tb -> ConditionalExpression[-Pi + ArcTan[1/Sqrt[35]] +
2*Pi*C[1], Element[C[1], Integers]]},
{x -> Sqrt[35]/2, tb -> ConditionalExpression[-ArcTan[1/Sqrt[35]] +
2*Pi*C[1], Element[C[1], Integers]]}}


but that's hardly surprising...).

More interesting is the general case (i. e. non-constant Ra, Rb, ta) :

sage: var("Ra, Rb, ta");
sage: Sys = [Ra*cos(ta)+Rb*cos(tb)==x, Ra*sin(ta)+Rb*sin(tb)==0]


Maxima can't solve this (no surprise...) :

sage: solve(Sys, [x, tb])
[Ra*cos(ta) + Rb*cos(tb) == x, Ra*sin(ta) + Rb*sin(tb) == 0]
sage: solve(Sys, [x, tb], to_poly_solve="force")
[Ra*cos(ta) + Rb*cos(tb) == x, Ra*sin(ta) + Rb*sin(tb) == 0]


but Sympy (and Mathematica) can :

sage: solve(Sys, [x, tb], algorithm="sympy")
[{tb: pi + arcsin(Ra*sin(ta)/Rb),
x: Ra*cos(ta) - sqrt(-Ra^2*sin(ta)^2/Rb^2 + 1)*Rb},
{tb: -arcsin(Ra*sin(ta)/Rb),
x: Ra*cos(ta) + sqrt(-Ra^2*sin(ta)^2/Rb^2 + 1)*Rb}]
sage: mathematica.Solve(Sys, [x, tb])
{{x -> Ra*Cos[ta] - Sqrt[Rb^2 - Ra^2*Sin[ta]^2],
tb -> ConditionalExpression[ArcTan[-(Sqrt[Rb^2 - Ra^2*Sin[ta]^2]/Rb),
-((Ra*Sin[ta])/Rb)] + 2*Pi*C[1], Element[C[1], Integers]]},
{x -> Ra*Cos[ta] + Sqrt[Rb^2 - Ra^2*Sin[ta]^2],
tb -> ConditionalExpression[ArcTan[Sqrt[Rb^2 - Ra^2*Sin[ta]^2]/Rb,
-((Ra*Sin[ta])/Rb)] + 2*Pi*C[1], Element[C[1], Integers]]}}


EDIT : Coaxing Maxima's default solver to give a solution.

You have a system of two equations with three unknowns ; one of the has to become a parametter of the solution(s). Choosing ta as a parameter allows us to write :

sage: S1 = Sys[1].solve(tb) ; S1
[tb == -arcsin(Ra*sin(ta)/Rb)]
sage: Sol = S1+Sys[0].subs(S1).solve(x) ; Sol
[tb == -arcsin(Ra*sin(ta)/Rb),
x == Ra*cos(ta) + sqrt(-Ra^2*sin(ta)^2/Rb^2 + 1)*Rb]


The (almost) equivalence of this answer and those given by Sympy ant the Wolfram engine is left as an (interesting !) exercise to the reader...

HTH,

more

This is good stuf. All I needed to know is the solve parameter algorithm="sympy" I have searched and searched, but can find no exaustive documentation for solve. Please point me to a source which covers all the optional parameters. Thanks again for your excelent help.

( 2023-01-14 13:14:03 +0100 )edit

From

solve?

• "algorithm" - string (default: 'maxima'); to use SymPy's solvers set this to 'sympy'. Note that SymPy is always used for diophantine equations. Another choice is 'giac'.
( 2023-01-14 13:53:46 +0100 )edit

All I needed to know is the solve parameter algorithm="sympy"

I beg to differ : a proficient use of Sage greatly benefits of judiciously using the "bells and whistles" added to the "basic" definition of its functions, which were added for a reason (more often than not, a good one, rather than a developer's whim...).

It is probably impossible to learn all those : you should aim to learn the existence of such options, and (frequently) use online help as needed...

can find no exaustive documentation for solve

solve?, then, for the various available solvers, maxima_calculus.solve?, import sympy then sympy.solve?, mathematica.Solve?, giac.solve? ; for fricas, you'd better use Fricas' doc on the Web (fricas.solve? is somewhat uninformative :-)...).

HTH,

( 2023-01-14 14:23:03 +0100 )edit

Ra=0.5;
Rb=3.0;
ta=pi/2;
var('x tb');
sol = solve([Ra*cos(ta)+Rb*cos(tb)==x,
Ra*sin(ta)+Rb*sin(tb)==0.0],x,tb,
algorithm='sympy')
sol

[{tb: -0.167448079219689, x: 2.95803989154981},
{tb: 3.30904073280948, x: -2.95803989154981}]

more