- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: ExSanQuoteOriginally posted by: Cuchulainne = [2,1,2,1,1,4,1,1,6,1,...,1,2n,1,...] wolf alpha

http://www.datasimfinancial.com

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

A follow-on is to solve Kepler's equation to 15 digits accuracy:M = E - eps * sin(E)with M = 0.8, eps = 0.2. Solve for E.How many ways? //E = 0.964333887695228

Last edited by Cuchulainn on January 19th, 2015, 11:00 pm, edited 1 time in total.

http://www.datasimfinancial.com

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: outrunand extra bonus points for eps=0.7Depends. What is the answer? Did you actually run it.

Last edited by Cuchulainn on January 28th, 2015, 11:00 pm, edited 1 time in total.

http://www.datasimfinancial.com

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

Using Newton-Raphson in this case is a disaster waiting to happen (for all the known reasons!!). It's probably the wrong method.For elliptic orbits (0 < eps < 1) you can see just by looking at the equation that the Banach fixed point method will work immediately And there is an analytical solution using series of Bessel functions of the first kind.I tested using fixed point and Differential Evolution with Least Squares (why not, cool?) to geteps, value0.2, 0.9643..0.7, 1.49810.8, 1.59960.9, 1.69321.0, 1.7785Hyperbolic orbits M = eps sinh(E) - E are even easier because the iteration function is always a contraction. // Kepler was probably not aware of the work of Leonardo of Pisa (Fibonacci to you) who used fixed point in the 12th century (!) to solve polynomial equations.

Last edited by Cuchulainn on January 30th, 2015, 11:00 pm, edited 1 time in total.

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

The initial guess X0 = Pi/2 + (sqrt(1+eps*(2*(M+eps)-Pi)) -1)/eps should do for Newton (using approximation of order 2 in E = Pi/2). Not quite a task I want to do by paper and pencil.

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: AVtThe initial guess X0 = Pi/2 + (sqrt(1+eps*(2*(M+eps)-Pi)) -1)/eps should do for Newton (using approximation of order 2 in E = Pi/2). Not quite a task I want to do by paper and pencil.NR seems to have issues in the following ranges. eps [0.5, 0.9]. M [0.0. 0.1] //The solution using Bessel functions (using Math Toolkit) give the same results as FP and DE (the series converges 0 < eps < 1). For eps = 1 we can get 4-digit accuracy bv taking more terms in the series.

Last edited by Cuchulainn on January 30th, 2015, 11:00 pm, edited 1 time in total.

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

Either allow complex numerics or use the absolute value of X0.

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

QuoteEither allow complex numerics or use the absolute value of X0. This is too ad-hoc. I want to get rid of seeds altogether. Too much messing around and time-consuming (e.g. when time is money).Just by looking at the equation f(E) we can bracket [a,b] it such that f(a)*f(b) < 0 and then use Bisection method1. as a high-order method.2. to bracket the solution to a low tol, apply Bisection and use the result as seed to NR. Now I get the desired results.QED.I hate NR. // 6 solutions to date.

Last edited by Cuchulainn on January 30th, 2015, 11:00 pm, edited 1 time in total.

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

I do not mind ad hoc. There where fixed parameters, no?I agree to reject complex numerics for usual code. But use order=3 in E=0, 0 <= M < 1/3. Solving for the (first) real solution will do as guess, it should be X0 = rho/eps + 2*(eps-1)/rho where rho = ((3*M+(-12/eps^2*d)^(1/2))*eps^2)^(1/3), d= -1/12*eps*(9*M^2*eps-8*eps^3+24*eps^2-24*eps+8)Splitting at M=1/3 should let converge Newton in 3 - 5 steps. There is only 1 solution for 0 < eps < 1, as the function E-eps*sin(E)-M increases w.r.t. E

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

QuoteOriginally posted by: AVtI do not mind ad hoc. There where fixed parameters, no?I agree to reject complex numerics for usual code. But use order=3 in E=0, 0 <= M < 1/3. Solving for the (first) real solution will do as guess, it should be X0 = rho/eps + 2*(eps-1)/rho where rho = ((3*M+(-12/eps^2*d)^(1/2))*eps^2)^(1/3), d= -1/12*eps*(9*M^2*eps-8*eps^3+24*eps^2-24*eps+8)Splitting at M=1/3 should let converge Newton in 3 - 5 steps. There is only 1 solution for 0 < eps < 1, as the function E-eps*sin(E)-M increases w.r.t. EI bet you did not spend Saturday evening working out this formulas by hand (place smile emoticon here) Did you generate it in Maple or something? And the rationale eludes me.This kind of approach only convinces me more that I don't like NR. I had something similar on a fixed-price fixed-income problem (I take on the project risk) to discover that the NR seed was word. As in this thread, I first bracket using Bisection method and then apply NR. It works really well. It's like interval arithmetic, the kind of stuff they put into the Apollo. ====When eps == 0 the seed explodes.

Last edited by Cuchulainn on February 1st, 2015, 11:00 pm, edited 1 time in total.

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

We can mitigate the dependence of NR on the choice of seeds by the bespoke (on this thread) Davidenko ODE (solve with Boost odeint) and Discrete Imbedding methods. The problematic cases M = 0.1, eps >= 0.8 are thus resolved.So, 8 solutions to date. 1. NR with seed2. Bisection standalone3. NR with Bisection pre-bracketer4. Fixed point5. Differential Evoulation + least squares6. Analytical solution (Bessel functions)7. Homotopy ..Davidenko ODE8. Discrete imbedding.Results 2-8 consistent with each other.

Last edited by Cuchulainn on January 31st, 2015, 11:00 pm, edited 1 time in total.

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

- Cuchulainn
**Posts:**61129**Joined:****Location:**Amsterdam-
**Contact:**

9. Brent minimization least squaresM, eps, E, error0.1, 0.7, 0.3205, O(1.0e-19)0.1, 0.8, 0.4427, O(1.0e-21)

Last edited by Cuchulainn on February 1st, 2015, 11:00 pm, edited 1 time in total.

http://www.datasim.nl

Every Time We Teach a Child Something, We Keep Him from Inventing It Himself

Jean Piaget

Yes, we have different attitudes for that, here I would prefer to havea good initial instead of a generic approach or naming methods. And soavoiding probably expansive methods.Since generations already considered that problem (and the best I haveread about is '2 steps') here is the one I like with Newton's method.For that I borrow a series solution for eps=1 from D. Cantrell to getestimates for the ugly case eps ~ 1 (try it for 1 - small, M ~ 0.25or close to 0).Note that by periodics one can reduce to M <= 2*Pi and moreover to thecase of M <= Pi with switching sign for eps by symmetry of sinus, butI only care for 0 <= eps <= 1.The code is generated by Maple, so clean it up in case. It is obvioushow to return the actual result=solution instead of number of steps.PS: edited the code region to provide 18 decimals instead of only 15.

Last edited by AVt on February 2nd, 2015, 11:00 pm, edited 1 time in total.

GZIP: On