Ask Your Question
0

How do I solve vector loop equations

asked 2023-01-13 21:38:11 +0200

updated 2023-01-14 14:27:55 +0200

Emmanuel Charpentier gravatar image

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 flag offensive close merge delete

Comments

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

brad.can@ntlworld.com gravatar imagebrad.can@ntlworld.com ( 2023-01-13 21:39:47 +0200 )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
tolga gravatar imagetolga ( 2023-01-14 06:49:36 +0200 )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.

brad.can@ntlworld.com gravatar imagebrad.can@ntlworld.com ( 2023-01-14 11:24:26 +0200 )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...

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-01-14 11:51:04 +0200 )edit

2 Answers

Sort by ยป oldest newest most voted
0

answered 2023-01-14 11:33:19 +0200

Emmanuel Charpentier gravatar image

updated 2023-01-14 14:43:11 +0200

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,

edit flag offensive delete link more

Comments

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.

brad.can@ntlworld.com gravatar imagebrad.can@ntlworld.com ( 2023-01-14 13:14:03 +0200 )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'.
achrzesz gravatar imageachrzesz ( 2023-01-14 13:53:46 +0200 )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,

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2023-01-14 14:23:03 +0200 )edit
0

answered 2023-01-14 12:11:22 +0200

achrzesz gravatar image

Abbreviated @Emmanuel answer

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}]
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2023-01-13 21:27:09 +0200

Seen: 125 times

Last updated: Jan 14 '23