Runge-Kutta-Verner method with error-estimate « Next Oldest | Next Newest »

 ▼ Jean-Michel Member Posts: 191 Threads: 41 Joined: Jun 2007 05-30-2007, 04:09 PM ```Hi! I've found on the web the coefficients of a Runge-Kutta-Verner method of order 8 embedded within 9th-order. Here is a program that uses these formulae for the HP-48: ------------------------- 'RKV8E9' ---------------------------- << OVER 1 SWAP START DUP 5 PICK ->NUM 4 PICK * DUP 12 / 3 PICK + 6 PICK ->NUM 5 PICK * 2 * OVER + 27 / 3 PICK + 6 PICK ->NUM 5 PICK * DUP2 3 * + 24 / 4 PICK + 7 PICK ->NUM 6 PICK * DUP 2.23331697733 * ROT -2.39805710715 * + 3 PICK .624672095524 * + 4 PICK + 7 PICK ->NUM 6 PICK * DUP .245675793931 * 3 PICK .273953453873 * + 4 PICK 4.36700683814E-2 * + 5 PICK + 8 PICK ->NUM 7 PICK * DUP 7.02490360993E-3 * ROT 1.34776896111E-2 * - ROT .181531822412 * + 3 PICK 6.16216474034E-2 * + 4 PICK + 7 PICK ->NUM 6 PICK * DUP .341657217459 * 3 PICK .250935375134 * + 4 PICK 7.40740740741E-2 * + 5 PICK + 8 PICK ->NUM 7 PICK * DUP -.03515625 * 3 PICK .340504422039 * + 4 PICK .120433077961 * + 5 PICK .07421875 * + 6 PICK + 9 PICK ->NUM 8 PICK * DUP -.296296296296 * 3 PICK 16 / - 4 PICK .310705427943 * + 5 PICK .305035312798 * + 6 PICK 7.63888888889E-2 * + 7 PICK + 10 PICK ->NUM 9 PICK * DUP -.26063186736 * 3 PICK 7.27205419732E-2 * + 4 PICK .011746330035 * - 5 PICK .37852828889 * + 7 PICK 7.11293665317E-2 * + 8 PICK + 11 PICK ->NUM 10 PICK * DUP -328.691358025 * 3 PICK 241.543474144 * - 4 PICK 626.941484898 * + 5 PICK 113.719201869 * + 6 PICK 413.485511011 * + 7 PICK 574.436392562 * - 8 PICK 8.14163971388 * - 9 PICK + 12 PICK ->NUM 11 PICK * DUP 1.80288461538E-3 * 3 PICK 2.49192782526 * + 4 PICK 7.69118880716E-2 * - 5 PICK .690428248367 * - 6 PICK .228863388685 * + 7 PICK 1.9030978898 * - 8 PICK .69337350173 * + 9 PICK .087803759282 * + 10 PICK + 13 PICK ->NUM 12 PICK * DUP 1.77777777778 * 3 PICK 2.35042735043E-2 * - 4 PICK 4.86229819563 * - 5 PICK 4.88888888889 * - 6 PICK 2.27160493827 * + 7 PICK 6.22916666667 * - 8 PICK 7.48550699458 * + 9 PICK 5.5746781906 * + 10 PICK .105709876543 * - 11 PICK + 14 PICK ->NUM 13 PICK * OVER 26.25 / 4 PICK 1920 / - 5 PICK .564373897707 * + 6 PICK .867950955719 * - 7 PICK .701048612636 * + 8 PICK .125460177737 * - 9 PICK 7.51961665302E-2 * + 10 PICK .272718881509 * - 11 PICK 5.46035999955E-2 * + 12 PICK + 15 PICK ->NUM 14 PICK * DUP 5.40772532189 * 4 PICK .422317596567 * + 5 PICK 2.87223506108E-3 * - 6 PICK 18.706283702 * - 7 PICK 7.80529021378 * + 8 PICK 3.73079546162 * + 9 PICK 1.38102723196 * + 10 ROLL 6.81524922033 * + 10 ROLL 5.45684741856 * - 10 PICK .396401690533 * - 11 PICK + 14 PICK ->NUM 13 PICK * -5.54761904762E-2 * SWAP .36 * - OVER 3.21428571429E-2 * + 3 PICK .12 * + 4 PICK 5.76923076923E-4 * + 5 PICK 1.05025641026 * + 6 PICK 1.05 * - 7 PICK .56 * + 8 PICK .315 * - 9 PICK .0175 * + 1 0 PUT 'ERROR' IFERR STO+ THEN STO END 3.21428571429E-2 * SWAP .342857142857 * + SWAP 4.12087912088E-4 * + SWAP .750183150183 * + SWAP .717857142857 * - SWAP .72380952381 * + SWAP .192857142857 * - SWAP 6.13095238095E-2 * + + NEXT >> ( 2368.5 bytes / #57545d ) *** PURGE 'ERROR' before the first execution. 4 inputs are needed: level 4: a program ( or its name ) that computes the derivative of the (n+1)-vector [ x y1 ... yn ] Since dx/dx = 1 the first component is always 1 level 3: the stepsize h level 2: the number of steps N level 1: the initial value [ x y1 ... yn ] -For example, to solve the system dy/dx = -yzu y(0)=1 dz/dx = x(y+z-u) z(0)=1 du/dx = xy-zu u(0)=2 with h = 0.1 , N = 10: << OBJ-> DROP -> X Y Z U << 1 Y Z U * * NEG Y Z + U - X * X Y * Z U * - 4 ->ARRY >> >> ENTER 0.1 ENTER 10 ENTER [ 0 1 1 2 ] and press the [RK8E9] key -It yields ( in 137 seconds ) in level 1 [ 1 0.258207906449 1.15762398084 0.842178311722 ] and 'ERROR' contains [ 0 2.7631E-11 3.27392E-11 7.636E-11 ] Regards, JMB. ``` ▼ Les Wright Posting Freak Posts: 1,368 Threads: 212 Joined: Dec 2006 06-02-2007, 10:32 PM JM, this is an excellent bit of work as usual. Your sample problem is somewhat faster on my 49G+, and I would expect similar performance on the 50G. I find in some examples I have worked with that the error estimate is a little optimistic. This is understandable--users will know that the error estimates of the Runge-Kutta methods are based on the higher derivative terms of Taylor series expansions, and the assumption that these derivative values remain roughly the same across each subinterval is often not a reasonable one. Moreover, the estimated error is typically used in more sophisticated routines to adjust stepsize--do better than a desired tolerance at the step, increase the stepsize, do worse, decrease it. Your routine raises some interesting questions about the pros and cons using higher order methods vs. repeated application of lower order methods. In which case is the computational work load best balanced by the performance? I believe that the 48 and 49 series calcs use RKF45 as the basis of an adaptive stepsize routine. One day, perhaps after all of the hubbub regarding the 35s dies down , I will post some thoughts on this, and perhaps issue a challenge on how to best find results to some challenging nonstiff IVPs. Les ▼ Jean-Michel Member Posts: 191 Threads: 41 Joined: Jun 2007 06-04-2007, 05:37 PM ```Hi, the error-estimates are algebraically added after each step. Adding the magnitudes could be a less optimistic alternative: Simply add OBJ-> EVAL ->LIST ABS OBJ-> ->ARRY just before 1 0 PUT 'ERROR' near the end I was worried when I saw that some coefficients of this Runge-Kutta-Verner method were of the order of 600: I thought it was going to produce great roundoff errors, but fortunately, the corresponding k-value ( namely k12 ) is weighted with a small 3/7280 and indeed, the results are usually satisfactory, even in the last decimal: For instance: y'= -2xy , y(0) = 1 h = 0.1 gives y(1) = 0.367879441185 h = 0.05 ----- y(1) = 0.367879441171 ( exact ) h = 0.025 ---- y(1) = 0.367879441171 ( still exact! ) At least, 'RKV8E9' gives an idea of the obtained accuracy with 16 evaluations of the function per step, instead of 11+2*11=33 evaluations if one uses 'RK8' with h and then with h/2. This may be an advantage if the functions are complicated. In "Numerical Recipes" they don't seem to like high-order Runge-Kutta formulas. On the other hand, they praise Bulirsch-Stoer methods ( "Runge-Kutta is for ploughing the fields, Burlisch-Stoer is a high-strung racehorse" ) and it is a little contradictory: though quite fantastic, these methods require more than 11 evaluations per step to achieve an 8th-order formula! Moreover, roundoff errors are also greater. I've always been fascinated by Runge-Kutta methods, especially if they are compared to the Taylor's method and the huge amount of calculus they require! And the size of the non-linear systems that must be solved to find high-order formulae! It seems like a touch of magic! Of course, the classical 'RK4' is probably enough for an HP-41, but we can try more accurate formulae with a 49G or 5OG. In fact, I wrote 'RK8' and 'RKV8E9' for completeness, perhaps it will be useful. However, I must say that, even now, I still prefer the HP-41 programming language with its direct arithmetic in the 4-level stack, and all the trickeries we can ( must? ) find to create neat programs. Speed may be a problem, but thanks to Warren Furlows and his excellent V41, HP-41 programs can run faster than HP-49 programs! But I stop here my "philosophical' point of view... Regards, JMB. ```

 Possibly Related Threads... Thread Author Replies Views Last Post A fast Bernoulli Number method for the HP Prime Namir 16 3,885 11-22-2013, 04:46 PM Last Post: Namir MS advert shows spreadsheet with obvious error BruceH 3 1,555 11-14-2013, 09:50 AM Last Post: Bill (Smithville, NJ) HP Prime: Rounding error in determinant Stephan Matthys 3 1,183 10-25-2013, 09:29 PM Last Post: Walter B Prime Error or Mine? toml_12953 12 2,555 10-22-2013, 10:35 AM Last Post: toml_12953 Explaination on How to Reset Caculator in Users guie: error Harold A Climer 5 1,770 10-15-2013, 02:11 AM Last Post: cyrille de Brébisson Repair of HP-34C - Error 0 and Error 9 Jeff Kearns 3 1,243 10-11-2013, 12:29 PM Last Post: Randy A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 5,055 09-13-2013, 01:58 PM Last Post: Pier Aiello Do You Think a Listing Error Was Made? Jim Johnson 13 2,832 09-04-2013, 09:23 AM Last Post: Eddie W. Shore HP Prime: Error reports : Here Patrice 111 15,204 07-24-2013, 05:52 PM Last Post: Thomas Klemm Downhill Simplex Method (Nelder Mead) for WP-34S Marcel Samek 0 658 07-07-2013, 03:05 AM Last Post: Marcel Samek

Forum Jump: