Kepler's equation (c1) 2*\a = alpha-epsilon*sin(alpha) (d1) 2 A = alpha - sin(alpha) epsilon was solved by Lagrange with an infinite series of Bessel functions (c6) alpha = 2*sum(bessel_j[n](epsilon*n)*sin(2*n*\a)/n,n,1,inf)+2*\a inf ==== bessel_j (epsilon n) sin(2 n A) \ n (d6) alpha = 2 > ------------------------------- + 2 A / n ==== n = 1 which is equivalent to the double sum alpha = 2*'sum((-1)^k*('sum(epsilon^(n+k)*(n-k)^(n+k-1)*sin(2*(n-k)*\a)/(2^(n+k)*n!),n,0,k))/k!,k,1,inf)+2*\a k-1 ==== n + k n + k - 1 k \ epsilon (n - k) sin(2 (n - k) A) (- 1) > ---------------------------------------------- inf / n + k ==== ==== 2 n! \ n = 0 alpha = 2 > ----------------------------------------------------------- / k! ==== k = 1 + 2 A This converges very poorly as |epsilon| nears 1. Although Lambert_w seems inadequate to solve Kepler's directly, we have alpha = 2*\a+im((w+1)*'sum(\p[k](w)*z^k/k!,k,0,inf)) inf k ==== P (w) z \ k alpha = 2 A + im((w + 1) > --------) / k! ==== k = 0 where w:=lambert_w(-(epsilon*%e^(-2*%i*\a)/2)) - 2 %i A epsilon %e w := lambert_w(- ------------------) 2 and z:=-epsilon^2/(w*(2*w+2)^2) 2 epsilon z = - ------------, 2 w (2 w + 2) w P (w) = -----, 0 w + 1 and P[>0](w) are (not obviously) polynomials \p[k]:=funmake('lambda,['([w]),factor(subst(lambert_w=lambda([[y]],w),exp(-2*k*w)*diff(exp(2*k*lambert_w(-epsilon*exp(2*%i*y)/2)),epsilon,k)/2/k/w^k*epsilon^k*(w+1)^(2*k-1)))]) P (w) = 1 1 P (w) = 3 w + 2 2 2 P (w) = 20 w + 26 w + 9 3 3 2 P (w) = 210 w + 404 w + 273 w + 64 4 4 3 2 P (w) = 3024 w + 7692 w + 7672 w + 3524 w + 625 5 5 4 3 2 P (w) = 55440 w + 175272 w + 230422 w + 156374 w + 54505 w + 7776 6 6 5 4 3 2 P (w) = 1235520 w + 4667952 w + 7605084 w + 6801904 w + 3507828 w 7 + 985830 w + 117649 ... Unfortunately, the radius of convergence is even worse, due to the k^(k-1) growth of the P[k] coefficients. But if we take exactly the same number of terms of twice the Lambert sum minus the Bessel sum, alpha = -'sum((-1)^k*('sum(epsilon^(n+k)*(n-k)^(n+k-1)*sin(2*(n-k)*\a)/(2^(n+k-1)*n!),n,0,k))/k!,k,1,kmax)+2*\a+2*im((w+1)*'sum(\p[k](w)*z^k/k!,k,0,kmax)) k-1 ==== n + k n + k - 1 k \ epsilon (n - k) sin(2 (n - k) A) (- 1) > ---------------------------------------------- kmax / n + k - 1 ==== ==== 2 n! \ n = 0 alpha = - > ----------------------------------------------------------- / k! ==== k = 1 kmax k ==== P (w) z \ k + 2 A + 2 im((w + 1) > --------) / k! ==== k = 0 the errors magically cancel, and we need only kmax ~ 13 terms for double precision, even when epsilon = 1, which, with alpha=y-x+pi/2, A = pi/4-x, lets us solve the descending tilted cosine wave y + x = cos(y - x) for y(x), an interesting involution. --rwg ----------- hbaker>FYI -- For those who 1) have a mortgage; that 2) isn't underwater; and 3) still have a job & can qualify for refinancing: "Optimal Mortgage Refinancing: A Closed Form Solution" http://www.economics.harvard.edu/faculty/laibson/files/Mortgage%20refinancin... http://www.nber.org/papers/w13487

An actual mortgage refinancing calculator based on the formulae from this paper: http://zwicke.nber.org/refinance/index.py