用 (wx)Maxima 求解方程:控制堆栈耗尽

Solving Equations with (wx)Maxima: Control stack exhausted

用 (wx)Maxima 求解方程:控制堆栈已耗尽

我正在尝试用 (wx)Maxima 求解方程:制定方程,然后让它插入变量并求解缺失变量的方程。但我很难过。最后一行有问题:

 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.

这是我的代码:

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);

如果我计算它 "directly",它仍然有效:

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

很遗憾,unit 软件包已经有一段时间没有更新了。我建议改用包 ezunits,其中量纲用反引号表示。要求解方程式,请尝试 dimensionally,它会通过一些旋转来帮助其他具有量纲的函数,例如dimensionally (solve (...))。 (请注意,dimensionally 未记录在案,对于不足之处,我深表歉意。)

我对你的程序做了一些修改,删除了一些不需要的东西,还使用了有理数而不是浮点数;与浮点数相比,Maxima 通常更适合有理数和整数。这是程序:

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 (%));

这是我得到的输出。请注意,我将解代回方程以验证它们。它看起来像预期的那样工作。

(%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

编辑。关于将 kg*m/s^2 转换为 N,您可以应用双反引号运算符。例如:

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

顺便说一下,要转换回浮点数,您可以应用 float:

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

FN 转换为 N 有点复杂,因为它是一个更复杂的表达式,尤其是因为 H 还没有附加单位。一些检查似乎显示 H 的单位必须是 kg*m/s^2。我将应用 declare_units 来说明 H 的单位是什么。然后我将 FN 转换为 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

符号qty(H)表示H的未指定数量。也可以只 subst(H = 100 ` kg*m/s^2, FN)(或任何数量,不只是 100)然后从那里开始。