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

user_4999991 0
Viewed 83

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 * 2 + ...

that approximately solves the

diff(A, zeta) = 0

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

polynomial-approximations wxmaxima maxima numeric
2 Answers

(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]) $
user_8710960

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]);
user_8710960
Related