CODEWITHSUNDEEP
maximawxmaxima
Dux
Dux

Reputation: 101

Solving Equations with (wx)Maxima: Control stack exhausted

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

Upvotes: 1

Views: 435

Answers (1)

Robert Dodier
Robert Dodier

Reputation: 17576

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.

Upvotes: 1

Related Questions