问题描述
对于我的小 FLOSS project，我想近似计算点接触的最大剪应力的 Green et al. 方程：
绘制时应该是这样的
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 + ...
大致解决了
diff(A,zeta) = 0
0 < ? < 0.5 和 0 < ? < 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,1);
(%o4) 0.4643131929806135
这里我将绘制不同 nu 值的近似解：
(%i5) plot2d (find_root (dAdzeta,1),[nu,0.5]) $
让我们将其与方程一起绘制。 10 这是格林在论文中推导出的近似值：
(%i6) plot2d ([find_root (dAdzeta,0.38167 + 0.33136*nu],0.5]) $
(2) 我研究了一些不同的方法来获得象征性的解决方案，这里有一些可能可行的方法。请注意，这也是一个近似值，因为它源自泰勒级数。你必须看看它是否足够好。
找到 acot 的低阶泰勒级数并将其插入 dAdzeta。
(%i7) acot_approx: taylor (acot(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,0.38167 + 0.33136*nu,rhs(zeta_max[3])],0.5]) $
这是一种不同的方法，它是通过 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,5);
(%o4) [0.0,0.1,0.2,0.3,0.4,0.5]
(%i5) makelist (find_root (dAdzeta,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],[1.094924180450318e-4])))]$
(%o9) done
不确定哪些部分最相关，所以我只返回了一些东西。可以通过 assoc 提取项目。这里我将提取构造函数。
(%i10) assoc ('f,foo);
X - 0.25
(%o10) lambda([X],block([Xtilde : ------------------,0.1707825127659933
2 3
X1 : [1,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,%o10],0.5]);