How to find a polynomial as an approximate solution to a nonlinear equation?

Question

For my small FLOSS project, I want to approximate the Green et al. equation for maximum shear stress for point contact:

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;

Now to find the maximum _max I differentiate the above equations against :

diff(A, zeta);

trying to solve the derivative for :

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

I ended up with a multipage equation that I can't actually use or test.

Now I was wondering if I can find the polynomial:

_max = a + b *  + c * ² + ...

that approximately solves the

diff(A, zeta) = 0

equation for 0 < < 0.5 and 0 < < 1.

score 1 · Answer 1 · answered Mar 11 '21 at 17:58

(1) Probably the first thing to try is just to solve diff(A, zeta) = 0 numerically (via find_root in this case). Here is an approximate solution for one value of 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

Here I'll plot the approximate solution for different values of nu:

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

Let's plot that together with Eq. 10 which is the approximation derived in the paper by Green:

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

(2) I looked at some different ways to get to a symbolic solution and here is something which is maybe workable. Note that this is also an approximation since it's derived from a Taylor series. You would have to look at whether it works well enough.

Find a low-order Taylor series for acot and plug it into 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

The approximate dAdzeta is a cubic polynomial in zeta, so we can solve it. The result is a big messy expression. The first two solutions are complex and the third is real, so I guess that's the one we want.

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

I tried some ideas to simplify the solution, but didn't find anything workable. Whether it's usable in its current form, I'll let you be the judge. Plotting the approximate solution along with the other two seems to show they're all pretty close together.

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

just as a small note, equation 10 in Green's paper is for maximum von Mises, not the max shear stress as far as I have understood. The rest of the instructions I'm following. — Foad S. Farimani, Mar 11 '21 at 22:30
BTW the solution doesn't have to be symbolic, actually. If Maxima had something like a list comprehension, as Python does, and decent polynomial curve fitting, this problem could be neatly solved. — Foad S. Farimani, Mar 11 '21 at 22:51
Foad, glad to hear it is helping you. About list comprehension, you can try either `makelist(f(x), x, mylist)` or `map(f, mylist)`. About polynomial fitting, I know there is a cubic spline `cspline`, and probably other polynomial approximation stuff elsewhere in maxima/share. It would be straightforward to put together a new package if needed. — Robert Dodier, Mar 11 '21 at 23:31
If we could have something like MATLAB's `polyfit()` that would be the ideal answer to my question I guess. — Foad S. Farimani, Mar 12 '21 at 16:43

score 1 · Answer 2 · answered Mar 13 '21 at 06:05

Here's a different approach, which is to calculate some approximate values by find_root and then assemble an approximation function which is a cubic polynomial. This makes use of a little function I wrote named polyfit. See: https://github.com/maxima-project-on-github/maxima-packages/tree/master/robert-dodier and then look in the polyfit folder.

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

Not sure which pieces are going to be most relevant, so I just returned several things. Items can be extracted via assoc. Here I'll extract the constructed function.

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

Plotting the function shows it is close to the values returned by find_root.

(%i12) plot2d ([find_root (dAdzeta, zeta, 0, 1), %o10], [nu, 0, 0.5]);

you are awesome. will look into this tonight. – Foad S. Farimani Mar 13 '21 at 06:22 — Foad S. Farimani, Mar 13 '21 at 06:22
Glad to help, it is fun to solve problems! – Robert Dodier Mar 13 '21 at 20:41 — Robert Dodier, Mar 13 '21 at 20:41

How to find a polynomial as an approximate solution to a nonlinear equation?

2 Answers2