numeric - 如何找到多项式作为非线性方程的近似解？

Question

给我的小FLOSS project ，我想近似Green et al.点接触的最大剪应力方程:

that should looks like this when plotted

the same equation in Maxima:

A: (3 / 2 / (1 + zeta^2) - 1 - nu + zeta * (1 + nu) * acot(zeta)) / 2;

现在求最大值 𝜏max我将上述方程与 𝜁 进行微分:

diff(A, zeta);

试图求解 𝜁 的导数:

solve(diff(A, zeta), zeta); 

我最终得到了一个我无法实际使用或测试的多页方程式。
现在我想知道是否可以找到多项式:

𝜁max = a + b * 𝜈 + c * 𝜈2 + ...

that approximately solves the

diff(A, zeta) = 0

0 < 𝜈 < 0.5 的方程和 0 < 𝜁 < 1 .

Answer 1

(1) 可能首先要尝试的就是解决diff(A, zeta) = 0数字(通过 find_root 在这种情况下)。这是 nu 的一个值的近似解:

(%i2) A: (3 / 2 / (1 + zeta^2) - 1 - nu + zeta * (1 + nu) * acot(zeta)) / 2;
                                         3
        (nu + 1) zeta acot(zeta) + ------------- - nu - 1
                                          2
                                   2 (zeta  + 1)
(%o2)   -------------------------------------------------
                                2
(%i3) dAdzeta: diff(A, zeta);
                             (nu + 1) zeta      3 zeta
       (nu + 1) acot(zeta) - ------------- - ------------
                                   2              2     2
                               zeta  + 1     (zeta  + 1)
(%o3)  --------------------------------------------------
                               2
(%i4) find_root (subst ('nu = 0.25, dAdzeta), zeta, 0, 1);
(%o4)                  0.4643131929806135

在这里，我将绘制 nu 不同值的近似解。 :

(%i5) plot2d (find_root (dAdzeta, zeta, 0, 1), [nu, 0, 0.5]) $

让我们将其与等式一起绘制。 10 这是格林在论文中得出的近似值:

(%i6) plot2d ([find_root (dAdzeta, zeta, 0, 1), 0.38167 + 0.33136*nu], [nu, 0, 0.5]) $

(2) 我研究了一些不同的方法来获得象征性的解决方案，这里有一些可能可行的方法。请注意，这也是一个近似值，因为它源自泰勒级数。你必须看看它是否足够好。
找到 acot 的低阶泰勒级数并将其插入 dAdzeta .

(%i7) acot_approx: taylor (acot(zeta), zeta, 1/2, 3);
                             1              1 2              1 3
                   4 (zeta - -)   8 (zeta - -)    16 (zeta - -)
                             2              2                2
(%o7)/T/ atan(2) - ------------ + ------------- + -------------- + . . .
                        5              25              375
                                                         
(%i8) dAdzeta_approx: subst (acot(zeta) = acot_approx, dAdzeta);
         (25 atan(2) - 10) nu + 25 atan(2) - 34
(%o8)/T/ --------------------------------------
                           50
                         1                            1 2
   (80 nu + 104) (zeta - -)   (320 nu + 1184) (zeta - -)
                         2                            2
 - ------------------------ + ---------------------------
             125                          625
                            1 3
   (640 nu + 11584) (zeta - -)
                            2
 - ---------------------------- + . . .
               9375

近似 dAdzeta是 zeta 中的三次多项式，所以我们可以解决它。结果是一个大杂乱的表情。前两个解决方案很复杂，第三个是真实的，所以我想这就是我们想要的。

(%i9) zeta_max: solve (dAdzeta_approx = 0, zeta);
<large mess omitted here>

(%i10) grind (zeta_max[3]);
zeta = ((625*sqrt((22500*atan(2)^2+30000*atan(2)-41200)*nu^4
                   +(859500*atan(2)^2-1878000*atan(2)+926000)
                    *nu^3
                   +(9022725*atan(2)^2-15859620*atan(2)+7283316)
                    *nu^2
                   +(15556950*atan(2)^2-36812760*atan(2)
                                       +19709144)
                    *nu+7371225*atan(2)^2-22861140*atan(2)
                   +17716484))
     /(256*(10*nu+181)^2)
     +((3*((9375*nu+9375)*atan(2)+4810*nu+6826))/(1280*nu+23168)
      -((90*nu+549)*(1410*nu+4281))/((10*nu+181)*(80*nu+1448)))
      /6+(90*nu+549)^3/(27*(10*nu+181)^3))
     ^(1/3)
     -((1410*nu+4281)/(3*(80*nu+1448))
      +((-1)*(90*nu+549)^2)/(9*(10*nu+181)^2))
      /((625*sqrt((22500*atan(2)^2+30000*atan(2)-41200)*nu^4
                   +(859500*atan(2)^2-1878000*atan(2)+926000)
                    *nu^3
                   +(9022725*atan(2)^2-15859620*atan(2)+7283316)
                    *nu^2
                   +(15556950*atan(2)^2-36812760*atan(2)
                                       +19709144)
                    *nu+7371225*atan(2)^2-22861140*atan(2)
                   +17716484))
       /(256*(10*nu+181)^2)
       +((3*((9375*nu+9375)*atan(2)+4810*nu+6826))
        /(1280*nu+23168)
        -((90*nu+549)*(1410*nu+4281))
         /((10*nu+181)*(80*nu+1448)))
        /6+(90*nu+549)^3/(27*(10*nu+181)^3))
       ^(1/3)+(90*nu+549)/(3*(10*nu+181))$

我尝试了一些想法来简化解决方案，但没有找到任何可行的方法。它是否可以以目前的形式使用，我会让你做判断。将近似解与其他两个一起绘制似乎表明它们都非常接近。

(%i18) plot2d ([find_root (dAdzeta, zeta, 0, 1),
                0.38167 + 0.33136*nu,
                rhs(zeta_max[3])], 
               [nu, 0, 0.5]) $