对于我的小型FLOSS 项目，我想近似Green 等人。点接触的最大剪应力方程：
绘制时应该是这样的
Maxima 中的相同方程：
A: (3 / 2 / (1 + zeta^2) - 1 - nu + zeta * (1 + nu) * acot(zeta)) / 2;
现在要找到最大值max，我将上述方程与 区分：
diff(A, zeta);
试图求解 的导数：
solve(diff(A, zeta), zeta);
我最终得到了一个我无法实际使用或测试的多页方程。
现在我想知道是否可以找到多项式：
最大值= a + b * + c * 2 + ...
这大约解决了
diff(A, zeta) = 0
0 < < 0.5和的方程0 < < 1。
这是一种不同的方法，它是通过计算一些近似值find_root，然后组装一个近似函数，它是一个三次多项式。这利用了我写的一个小函数，命名为polyfit. 请参阅：https ://github.com/maxima-project-on-github/maxima-packages/tree/master/robert-dodier然后查看polyfit文件夹。
(%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) nn: makelist (k/10.0, k, 0, 5);
(%o4) [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]
(%i5) makelist (find_root (dAdzeta, zeta, 0, 1), nu, nn);
(%o5) [0.3819362006941755, 0.4148794361988409,
0.4478096487716516, 0.4808644852928955, 0.5141748609122403,
0.5478684611102143]
(%i7) load ("polyfit.mac");
(%o7) polyfit.mac
(%i8) foo: polyfit (nn, %o5, 3) $
(%i9) grind (foo);
[beta = matrix([0.4643142407230925],[0.05644202066198245],
[2.746081069103333e-4],[1.094924180450318e-4]),
Yhat = matrix([0.3819365703555216],[0.4148782994206623],
[0.4478104992708994],[0.4808650578507559],
[0.5141738631047557],[0.5478688029774219]),
residuals = matrix([-3.696613460890674e-7],
[1.136778178534303e-6],
[-8.504992477509354e-7],
[-5.725578604010018e-7],
[9.97807484637292e-7],
[-3.418672076538343e-7]),
mse = 5.987630959972099e-13,Xmean = 0.25,
Xsd = 0.1707825127659933,
f = lambda([X],
block([Xtilde:(X-0.25)/0.1707825127659933,X1],
X1:[1,Xtilde,Xtilde^2,Xtilde^3],
X1 . matrix([0.4643142407230925],
[0.05644202066198245],
[2.746081069103333e-4],
[1.094924180450318e-4])))]$
(%o9) done
不确定哪些部分最相关，所以我只返回了几件东西。可以通过 提取项目assoc。这里我将提取构造函数。
(%i10) assoc ('f, foo);
X - 0.25
(%o10) lambda([X], block([Xtilde : ------------------, X1],
0.1707825127659933
2 3
X1 : [1, Xtilde, Xtilde , Xtilde ],
[ 0.4643142407230925 ]
[ ]
[ 0.05644202066198245 ]
X1 . [ ]))
[ 2.746081069103333e-4 ]
[ ]
[ 1.094924180450318e-4 ]
(%i11) %o10(0.25);
(%o11) 0.4643142407230925
绘制函数显示它接近返回的值find_root。
(%i12) plot2d ([find_root (dAdzeta, zeta, 0, 1), %o10], [nu, 0, 0.5]);
(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 这是 Green 在论文中得出的近似值：
(%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]) $