1

Solving Equations with (wx)Maxima: Control stack exhausted

I'm trying to solve equations with (wx)Maxima: formulate the equation, then let it insert the variables and solve the equation for the missing variable. But I'm having a hard time. Somehow it's having problems in the last line:

 Control stack exhausted (no more space for function call frames).
This is probably due to heavily nested or infinitely recursive function
calls, or a tail call that SBCL cannot or has not optimized away.

That's my code:

kill(all); 
load(physical_constants); 
load(unit); 
setunits([kg,m,s,N]); 
showtime: false;

α: 30*%pi/180; 
/*α: 30*°;*/ 
masse: 1000*kg; 
g: 9.80665*m/(s*s); 
b: 0.3*m; 
B: 0.5*m; 
L: 0.1*m;

F_g: masse*g; 
F_H: masse * g;

kill(S, x); 
S: solve(0=F_H-2*x*sin(α), x); 
S: assoc(x, S);

kill(H, x);
H: solve(0=-F_g+2*x, x);
H: assoc(x, H);

kill(Ly, x); 
Ly: solve(tan(α)=x/(B/2), x); 
Ly: assoc(x, Ly);

kill(FN, x); 
FN: solve(0=H*B/2-x*(L+Ly)+S*sin(α)*B/2+S*cos(α)*Ly, x); 
FN: assoc(x, FN);

If I calculate it "directly", it works though:

kill(all); 
load(physical_constants); 
load(unit); 
setunits([kg,m,s,N]); 
showtime: false;

kill(FN, x); 
FN: solve([α=30*%pi/180, H=196133/40*N,
           B=0.5*m, L=0.1*m, 
           Ly=sqrt(3)/12*m, S=196133/20*N,
           0=H*B/2-x*(L+Ly)+S*sin(α)*B/2+S*cos(α)*Ly], 
          [x, α, H,  B, L, Ly, S]); 
FN: assoc(x, FN[1]); 
FN: float(FN);

(FN)    1934473685/128529*N
Dux
  • 101
  • 4

1 Answers1

1

Unfortunately the unit package has not been updated in some time. I'll suggest to use instead the package ezunits, in which dimensional quantities are represented with a back quote. To solve equations, try dimensionally which goes through some gyrations to help other functions with dimensional quantities, e.g. dimensionally (solve (...)). (Note that dimensionally isn't documented, I'm sorry for the shortcoming.)

I've modified your program a little to remove some unneeded stuff and also to use rational numbers instead of floats; Maxima is generally more comfortable with rationals and integers than with floats. Here is the program:

linel: 65 $
load(ezunits) $

α: 30*%pi/180; 
masse: 1000`kg; 
g: rationalize(9.80665)`m/(s*s); 
b: 3/10`m; 
B: 5/10`m; 
L: 1/10`m;

F_g: masse*g; 
F_H: masse * g;

S: dimensionally (solve(0=F_H-2*x*sin(α), x)); 
S: assoc(x, S);

Ly: dimensionally (solve(tan(α)=x/(B/2), x));
Ly: assoc(x, Ly);

FN: dimensionally (solve(0=H*B/2-x*(L+Ly)+S*sin(α)*B/2+S*cos(α)*Ly, x));
FN: assoc(x, FN);

subst (x = S, F_H-2*x*sin(α));
subst (x = Ly, tan(α)=x/(B/2));
subst (x = FN, H*B/2-x*(L+Ly)+S*sin(α)*B/2+S*cos(α)*Ly);
ratsimp (expand (%));

and here is the output I get. Note that I substituted the solutions back into the equations to verify them. It looks like it worked as expected.

(%i2) linel:65
(%i3) load(ezunits)
(%i4) α:(30*%pi)/180
                               %pi
(%o4)                          ---
                                6
(%i5) masse:1000 ` kg
(%o5)                       1000 ` kg
(%i6) g:rationalize(9.80665) ` m/(s*s)
                      5520653160719109   m
(%o6)                 ---------------- ` --
                      562949953421312     2
                                         s
(%i7) b:3/10 ` m
                             3
(%o7)                        -- ` m
                             10
(%i8) B:5/10 ` m
                              1
(%o8)                         - ` m
                              2
(%i9) L:1/10 ` m
                             1
(%o9)                        -- ` m
                             10
(%i10) F_g:masse*g
                    690081645089888625   kg m
(%o10)              ------------------ ` ----
                      70368744177664       2
                                          s
(%i11) F_H:masse*g
                    690081645089888625   kg m
(%o11)              ------------------ ` ----
                      70368744177664       2
                                          s
(%i12) S:dimensionally(solve(0 = F_H-2*x*sin(α),x))
                      690081645089888625   kg m
(%o12)           [x = ------------------ ` ----]
                        70368744177664       2
                                            s
(%i13) S:assoc(x,S)
                    690081645089888625   kg m
(%o13)              ------------------ ` ----
                      70368744177664       2
                                          s
(%i14) Ly:dimensionally(solve(tan(α) = x/(B/2),x))
                                1
(%o14)                 [x = --------- ` m]
                            4 sqrt(3)
(%i15) Ly:assoc(x,Ly)
                              1
(%o15)                    --------- ` m
                          4 sqrt(3)
(%i16) FN:dimensionally(solve(0 = (H*B)/2-x*(L+Ly)
                                         +(S*sin(α)*B)/2
                                         +S*cos(α)*Ly,x))
                                 1                       1
(%o16) [x = (----------------------------------------- ` --)
             140737488355328 sqrt(3) + 351843720888320    2
                                                         s
                               2
 (351843720888320 sqrt(3) H ` s
                        3/2
 + 1150136075149814375 3    ` kg m)]
(%i17) FN:assoc(x,FN)
                            1                       1
(%o17) (----------------------------------------- ` --)
        140737488355328 sqrt(3) + 351843720888320    2
                                                    s
                               2
 (351843720888320 sqrt(3) H ` s
                        3/2
 + 1150136075149814375 3    ` kg m)
(%i18) subst(x = S,F_H-2*x*sin(α))
                                kg m
(%o18)                      0 ` ----
                                  2
                                 s
(%i19) subst(x = Ly,tan(α) = x/(B/2))
                           1         1
(%o19)                  ------- = -------
                        sqrt(3)   sqrt(3)
(%i20) subst(x = FN,(H*B)/2-x*(L+Ly)+(S*sin(α)*B)/2+S*cos(α)*Ly)
                           1        1
                    (- ---------) - --
                       4 sqrt(3)    10               1
(%o20) ((----------------------------------------- ` --)
         140737488355328 sqrt(3) + 351843720888320    2
                                                     s
                               2
 (351843720888320 sqrt(3) H ` s
                        3/2           H
 + 1150136075149814375 3    ` kg m) + -) ` m
                                      4
                            2
   690081645089888625   kg m
 + ------------------ ` -----
    281474976710656       2
                         s
(%i21) ratsimp(expand(%))
                                    2
                                kg m
(%o21)                      0 ` -----
                                  2
                                 s

EDIT. About converting kg*m/s^2 to N, you can apply the double back quote operator. For example:

(%i25) F_g `` N
                     690081645089888625
(%o25)               ------------------ ` N
                       70368744177664

By the way, to convert back to floats, you can apply float:

(%i26) float(%)
(%o26)                9806.649999999998 ` N

Converting FN to N is a little more involved, since it's a more complex expression, especially because of H which doesn't have units attached to it yet. Some inspection seems to show the units of H must be kg*m/s^2. I'll apply declare_units to say that's what are the units of H. Then I'll convert FN to N.

(%i27) declare_units(H,(kg*m)/s^2)
                              kg m
(%o27)                        ----
                                2
                               s
(%i28) FN `` N
             351843720888320 sqrt(3) qty(H)
(%o28) (-----------------------------------------
        140737488355328 sqrt(3) + 351843720888320
                                                3/2
                           1150136075149814375 3
                 + -----------------------------------------) ` N
                   140737488355328 sqrt(3) + 351843720888320
(%i29) float(%)
(%o29) (1.023174629940149 qty(H) + 10033.91548470256) ` N

The notation qty(H) represents the unspecified quantity of H. One could also just subst(H = 100 ` kg*m/s^2, FN) (or any quantity, not just 100) and go from there.

Robert Dodier
  • 16,905
  • 2
  • 31
  • 48
  • Thanks, load(ezunits) with dimensionally() does the trick. Though: Using ```*``` for the units instead of ``` ` ``` is easier on my keyboard. Is there a way to change the key? Also load(unit): ```Evaluation took 0.6240 seconds (0.6250 elapsed) using 219.953 MB.``` load(ezunits): ```Evaluation took 3.7280 seconds (3.7360 elapsed) using 1168.046 MB.``` And how do I get that setunits([kg,m,s,N]); behaviour? Could be in Newtons: ```(S) 196133/20 ` (kg*m)/s^2``` – Dux May 01 '20 at 09:19
  • Yes maxima likes fractions, but: ``` g: rationalize(9.80665)`m/s^2; g=5520653160719109/562949953421312=9.8066499999999994 grat: rat(9.80665)`m/s^2; grat=196133/20000=9.80665 ``` So what's the difference between rat() and rationalize() and why should I use the latter when it's just has bigger fractions which are more imprecise? Also, can you use g as a variable, since it also means gramm, or is that fixed with the ``` ` ``` for the units instead of ```*```? – Dux May 01 '20 at 09:19
  • Hi, thanks for your interest. I'll try to answer these one at a time. (1) About back quote versus asterisk: no, units must be denoted with back quote. That's necessary to distinguish nondimensional factors from units. (2) Yeah, loading ezunits is relatively slow. I'm thinking about how to make it faster. (I wrote the package.) (3) About converting to N, apply the double back quote. I've edited my answer to show that. – Robert Dodier May 01 '20 at 19:02
  • (4) About fractions, `rationalize` is the exact rational equivalent of a float. I guess that's probably too much precision for many cases. `ratsimp` and `rat` will try to find a simpler approximation. Note that `rat` creates a special representation of polynomials with rational coefficients; for just approximating floats, apply `ratsimp` instead. (5) ezunits doesn't restrict names for variables -- it won't complain if `g` is a variable in your program, but it is also a unit at the same time. This can lead to "gotcha" problems -- e.g. ``g: 98/10 ` m/s^2`` and then ``m1: 200 ` g`` makes a mess. – Robert Dodier May 01 '20 at 19:06
  • More about (5). I've thought about how to avoid problems like that, but haven't reached any conclusions yet. I think probably a simple strategy to avoid most problems is to quote (not evaluate) any symbols appearing on the right-hand side of the back quote. – Robert Dodier May 01 '20 at 19:08
  • Sure I'm interested. A FOSS math program? Sign me up! My bad: somehow the calcuation of H was lost(fixed in the original post): ``` kill(H, x); H: solve(0=-F_g+2*x, x); H: assoc(x, H); ``` So it should be half of F_g, and the unit should be Newton also. ```F_g `` N``` just defines it as N, not regarding what units it had before? That's not what I need. – Dux May 01 '20 at 20:21
  • ```g: 98/10 ` m/s^2; m1: 200 ` g;``` I though that was fixed by using the backquote instead of the asteriks. Than it's just harder to type, but breaks just as easily as "unit". Shouldn't it know that after the backquote it must be gram, since the var has a number in it, which cannot be a unit in the backquote? float() didn't always work, but I think with ezunits and dimensionally() it works now. If you're the author, maybe you're interested in another calculation I did, to see where the problems come from. It just worked with "unit", but ezunits only worked with the dimesionally() trick. – Dux May 01 '20 at 20:21
  • Re: "```F_g `` N``` just defines it as N, not regarding what units it had before?" -- no, the double back quote is a conversion operator, it converts the units of `F_g` to `N` (assuming the units are the same dimensionally). Yes, I'm interested in any other problems you have worked on. Maybe at this point it's appropriate to open a new question for other problems. About quoting `g` on the right-hand side -- yeah, I've bumped into that myself. I don't see an easy work around ... let me think about that. – Robert Dodier May 01 '20 at 20:53
  • For the record, I've written up my answer in a neater form and posted it to a blog: https://maxima-solved.blogspot.com/2020/05/solving-equations-with-dimensional.html Maybe that's helpful in some way. – Robert Dodier May 03 '20 at 21:47
  • Shall I send you the scripts to your gmail? Or what would be the best way of communicating? – Dux May 08 '20 at 10:49
  • Thanks for your interest. I think the best way is to post a message to the Maxima mailing list. You can subscribe to the list on this page: https://sourceforge.net/projects/maxima/lists/maxima-discuss That way others can participate in the discussion, contributing ideas and perhaps learning something that is useful to them. – Robert Dodier May 08 '20 at 14:11