SERVING THE QUANTITATIVE FINANCE COMMUNITY

 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

October 31st, 2010, 8:19 pm

That's very interesting. The algo is described in section 3.2 here, sorry in Dutch. What is noteworthy is they use Gram-Schmidt to produce an orthormal basis before LLL algo.LLL in EnglishHad a quick look at Constructive Dedekind reals looks like computing a real in an approximate interval. This looks cool. And now an upper limit How to LLL to do this?
Last edited by Cuchulainn on October 30th, 2010, 10:00 pm, edited 1 time in total.
 
User avatar
AVt
Posts: 1074
Joined: December 29th, 2001, 8:23 pm

exp(5) = [$]e^5[$]

November 1st, 2010, 8:46 pm

@wileysw:That can be done in Maple as well, though it does not have the command 'RootApproximant'.One defines a 'base' (which is one for transcendent input, not otherwise in general) bypowers and seeks for integer relations (with small numerical errors, thus depending onthe computational precision used), for which a variant of LLL is used.ApproxAlgebraicRelation:= proc(z,deg:osint, x)# z = any constant, # deg = desired degree for an approximate algebraic relation, # x = indeterminate to be used in the polynomiallocal w,u;#Digits:= Digits+3; # make your choice, Digits:= 2*Digits is finew:=Vector(deg+1, 'k -> z^(k-1)'); # = 1, z, z^2, z^3 ... w:=convert(%, list); # syntax needed for the following commandu:=IntegerRelations:-LinearDependency(evalf(w));0=sum( u[ i]*x^(i-1),i=1..deg+1); end proc;Using 14 digits (i.e. ~ hardware double precision) that replies by0 = -4491-13581*x+4778*x^2 for input z=Pi and variable x, degree=2.That says: Pi 'almost' satisfies that quadric. Within that code one can increase the accuracy - that leads to better results, but somelarger coefficients. And if one wants the 'most simple' one just increase degree, untilthe error falls below a limit (so there is a kind of trade off in combining degree andaccuracy (leading to longer expressions), what one wants to call 'most simple' ...).But that 'simple' approach is unlikely to produce those nice & curious results with veryhigh accuracy given at the Mathworld.@Cuchulainn:That 'rationalize' is degree=1, so can be done by LLL. However that is not needed at all.Essentially it (depending on computational precision) an iterated divison with remainder,until the remainder is 'killed' by the 'system epsilon' (~ precision). Silently you useit every day while coding ... (though it may differ from IEEE for other reasons).
Last edited by AVt on October 31st, 2010, 11:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

November 7th, 2010, 9:56 am

QuoteThat 'rationalize' is degree=1, so can be done by LLL. However that is not needed at all.Essentially it (depending on computational precision) an iterated divison with remainder,until the remainder is 'killed' by the 'system epsilon' (~ precision). Silently you useit every day while coding ... (though it may differ from IEEE for other reasons).If I understand correctly, the iteration converges monotonically until machine accuracy cannot see the difference?It is a one-sided estimate?
 
User avatar
AVt
Posts: 1074
Joined: December 29th, 2001, 8:23 pm

exp(5) = [$]e^5[$]

November 7th, 2010, 4:37 pm

Hm ... it really is like converting a decimal input to a machine numberas the read commands do (which depend on the rounding rules, which maybe what you mean by 'one-sided'), I once went through that (do not haveit hand). Practical it is an enhanced Euclidian algorithm. It stops, assoon as all bits are used (or likewise by 'system epsilon').And a IEEE float is just a special fraction & rational number.
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

December 3rd, 2010, 9:08 am

Don't know if these solutions have already been given:1. Compoundingexp(x) = lim (1+x/n)^n as n -> INFINITY2. y = exp(x) solves the differential equationdy/dx = y, y(0) = 1which you solve by Euler's fd method. Maybe that's the way Euler did it at the time?? It's first-order accurate, so you can determine the accuracy based on the interview constraints (neurons work slow).
Last edited by Cuchulainn on December 2nd, 2010, 11:00 pm, edited 1 time in total.
 
User avatar
ThinkDifferent
Posts: 659
Joined: March 14th, 2007, 1:09 pm

exp(5) = [$]e^5[$]

December 3rd, 2010, 12:47 pm

QuoteOriginally posted by: slopse=2.7 Ibsen Ibsen. (Ibsen was born in 1828, e = 2.718281828...) Ibsen? LOL The way I learned e = 2.718281828 in high school was that Lev Tolstoy was born in 1828.
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

January 28th, 2011, 1:26 pm

QuoteYou can get e^5~148 with sub-1% accuracy without pen and paper if you remember that ln(2)~0.693. Try it.Here's another way and you only have to remember that exp(0) = 11. Compute e == exp(1) using the fact that the exponential satisfies du/dx = u, u(0) = 1. Solve using Euler's method U(n+1) = (1+h)*U(n). Compute e ~ U(N) = (1+h)^N.Perhaps surprisingly, since this work on logarithms had come so close to recognising the number e, when e is first "discovered" it is not through the notion of logarithm at all but rather through a study of compound interest. In 1683 Jacob Bernoulli looked at the problem of compound interest and, in examining continuous compound interest, he tried to find the limit of (1 + 1/n)n as n tends to infinity. He used the binomial theorem to show that the limit had to lie between 2 and 3 so we could consider this to be the first approximation found to e. Also if we accept this as a definition of e, it is the first time that a number was defined by a limiting process. He certainly did not recognise any connection between his work and that on logarithms. 2. Now e^5 = (e^2) * (e^2) * e.
Last edited by Cuchulainn on January 29th, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
ExSan
Posts: 4554
Joined: April 12th, 2003, 10:40 am

exp(5) = [$]e^5[$]

February 3rd, 2011, 2:19 pm

QuoteOriginally posted by: CuchulainnHere's another way and you only have to remember that exp(0) = 11. Compute e == exp(1) using the fact that the exponential satisfies du/dx = u, u(0) = 1. Solve using Euler's method U(n+1) = (1+h)*U(n). Compute e ~ U(N) = (1+h)^N. Thanks to Cuchulainn for the hints ! Numerical Ordinary Differential EquationsFrom the Downloadable File extract program ExSan_EULER_ODE_EXP to your desktop and execute. (see updated post)
Last edited by ExSan on February 8th, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

February 4th, 2011, 8:23 am

A more accurate approx is to use the Cayley transform (aka Crank Nicolson)du/dx = u(U(n+1) - U(n))/h = (U(n+1) + U(n))/2or U(n+1) = alpha*U(n)alpha = (1 +h/2)/(1 - h/2)
Last edited by Cuchulainn on February 3rd, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
daveangel
Posts: 17031
Joined: October 20th, 2003, 4:05 pm

exp(5) = [$]e^5[$]

February 4th, 2011, 9:13 am

QuoteOriginally posted by: CuchulainnA more accurate approx is to use the Cayley transform (aka Crank Nicolson)du/dx = u(U(n+1) - U(n))/h = (U(n+1) + U(n))/2or U(n+1) = alpha*U(n)alpha = (1 +h/2)/(1 - h/2)neat
knowledge comes, wisdom lingers
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

February 7th, 2011, 1:34 pm

In general, it's better to solve for du/dx = -u, compute CN and then invert the answer! For the rest, it's the same as before.
Last edited by Cuchulainn on February 6th, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
ExSan
Posts: 4554
Joined: April 12th, 2003, 10:40 am

exp(5) = [$]e^5[$]

February 8th, 2011, 1:57 pm

IMPROVED EULER'S METHOD -Heun's Formula - Midpoint methodRegular Euler Method: exp( 5 ) ~ 129.1299 error = -12.99 % step = .1Improved Euler Method: exp( 5 ) ~ 147.2699 error = --0.77 % step = .1From the Downloadable File extract program ExSan_IMP_EULER_ODE_EXP to your desktop and execute. (see updated post)
Last edited by ExSan on February 8th, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

February 8th, 2011, 2:56 pm

QuoteOriginally posted by: CuchulainnIn general, it's better to solve for du/dx = -u, compute CN and then invert the answer! For the rest, it's the same as before.For systemsdU/dt = QU + rWe wish to compute a bounded solution as t -> infinity and this is determined by Q (See Varga 1962) and when the eigenvalues are < 0 then limit U(t) = Inverse(Q)*r as t -> infinity @exsanWhat is the order of accuracy h^n for your two schemes (1 and 2?) The convergence of CT to exp(-z) http://en.wikipedia.org/wiki/Pad%C3%A9_table
Last edited by Cuchulainn on February 7th, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
ExSan
Posts: 4554
Joined: April 12th, 2003, 10:40 am

exp(5) = [$]e^5[$]

February 9th, 2011, 1:40 pm

Runge - Kutta 4 OrderRK4 exp( 5 ) ~ 31.71365 error = -78.63 % step = 0.1RK4 exp( 5 ) ~ 28.37667 error = - -80.88 % step = 0.01Basic EULER Method is the most stable, so far, for ode du/dx = u , other methods are stable only for certain values of step (h)From the Downloadable File extract program ExSan_RK_EXP to your desktop and execute.
Attachments
ExSan_V.4.01.Y_RungeKutta_exp(x).zip
(1.63 MiB) Downloaded 46 times
Last edited by ExSan on February 8th, 2011, 11:00 pm, edited 1 time in total.
 
User avatar
Cuchulainn
Posts: 62395
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

exp(5) = [$]e^5[$]

March 5th, 2011, 12:35 pm

One scheme that might be useful are Exponential integrators esp. page 8.Have not tried this but I once did exponential fitting for this DE.
Last edited by Cuchulainn on March 4th, 2011, 11:00 pm, edited 1 time in total.
ABOUT WILMOTT

PW by JB

Wilmott.com has been "Serving the Quantitative Finance Community" since 2001. Continued...


Twitter LinkedIn Instagram

JOBS BOARD

JOBS BOARD

Looking for a quant job, risk, algo trading,...? Browse jobs here...


GZIP: On