1

使用 (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);

如果我“直接”计算它，它虽然有效：

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
4

1 回答 1

1

不幸的是，该unit软件包已经有一段时间没有更新了。我建议改用 package 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）然后从那里开始。

于 2020-04-30T03:58:23.317 回答