Terry Feagin's 10th order explicit Runge-Kutta method

Question

The following Julia code implements Terry Feagin's 10th order explicit Runge-Kutta method (a more accurate cousin of RK4). Though the structure of the code is quite simple (i.e. no cyclomatic complexity), it involves many high-precision coefficients and lengthy arithmetic expressions which bring its length to over 150 lines. The coefficients display no obvious pattern and cannot be computed at runtime.

With that said, are there any ways to simplify this function?

function rk108_step(f, x, y, h)
  const a0100 = BigFloat(1//10)

  const a0200 = BigFloat("-0.915176561375291440520015019275342154318951387664369720564660")
  const a0201 = BigFloat("1.45453440217827322805250021715664459117622483736537873607016")

  const a0300 = BigFloat("0.202259190301118170324681949205488413821477543637878380814562")
  const a0302 = BigFloat("0.606777570903354510974045847616465241464432630913635142443687")

  const a0400 = BigFloat("0.184024714708643575149100693471120664216774047979591417844635")
  const a0402 = BigFloat("0.197966831227192369068141770510388793370637287463360401555746")
  const a0403 = BigFloat("-0.0729547847313632629185146671595558023015011608914382961421311")

  const a0500 = BigFloat("0.0879007340206681337319777094132125475918886824944548534041378")
  const a0503 = BigFloat("0.410459702520260645318174895920453426088035325902848695210406")
  const a0504 = BigFloat("0.482713753678866489204726942976896106809132737721421333413261")

  const a0600 = BigFloat("0.0859700504902460302188480225945808401411132615636600222593880")
  const a0603 = BigFloat("0.330885963040722183948884057658753173648240154838402033448632")
  const a0604 = BigFloat("0.489662957309450192844507011135898201178015478433790097210790")
  const a0605 = BigFloat("-0.0731856375070850736789057580558988816340355615025188195854775")

  const a0700 = BigFloat("0.120930449125333720660378854927668953958938996999703678812621")
  const a0704 = BigFloat("0.260124675758295622809007617838335174368108756484693361887839")
  const a0705 = BigFloat("0.0325402621549091330158899334391231259332716675992700000776101")
  const a0706 = BigFloat("-0.0595780211817361001560122202563305121444953672762930724538856")

  const a0800 = BigFloat("0.110854379580391483508936171010218441909425780168656559807038")
  const a0805 = BigFloat("-0.0605761488255005587620924953655516875526344415354339234619466")
  const a0806 = BigFloat("0.321763705601778390100898799049878904081404368603077129251110")
  const a0807 = BigFloat("0.510485725608063031577759012285123416744672137031752354067590")

  const a0900 = BigFloat("0.112054414752879004829715002761802363003717611158172229329393")
  const a0905 = BigFloat("-0.144942775902865915672349828340980777181668499748506838876185")
  const a0906 = BigFloat("-0.333269719096256706589705211415746871709467423992115497968724")
  const a0907 = BigFloat("0.499269229556880061353316843969978567860276816592673201240332")
  const a0908 = BigFloat("0.509504608929686104236098690045386253986643232352989602185060")

  const a1000 = BigFloat("0.113976783964185986138004186736901163890724752541486831640341")
  const a1005 = BigFloat("-0.0768813364203356938586214289120895270821349023390922987406384")
  const a1006 = BigFloat("0.239527360324390649107711455271882373019741311201004119339563")
  const a1007 = BigFloat("0.397774662368094639047830462488952104564716416343454639902613")
  const a1008 = BigFloat("0.0107558956873607455550609147441477450257136782823280838547024")
  const a1009 = BigFloat("-0.327769124164018874147061087350233395378262992392394071906457")

  const a1100 = BigFloat("0.0798314528280196046351426864486400322758737630423413945356284")
  const a1105 = BigFloat("-0.0520329686800603076514949887612959068721311443881683526937298")
  const a1106 = BigFloat("-0.0576954146168548881732784355283433509066159287152968723021864")
  const a1107 = BigFloat("0.194781915712104164976306262147382871156142921354409364738090")
  const a1108 = BigFloat("0.145384923188325069727524825977071194859203467568236523866582")
  const a1109 = BigFloat("-0.0782942710351670777553986729725692447252077047239160551335016")
  const a1110 = BigFloat("-0.114503299361098912184303164290554670970133218405658122674674")

  const a1200 = BigFloat("0.985115610164857280120041500306517278413646677314195559520529")
  const a1203 = BigFloat("0.330885963040722183948884057658753173648240154838402033448632")
  const a1204 = BigFloat("0.489662957309450192844507011135898201178015478433790097210790")
  const a1205 = BigFloat("-1.37896486574843567582112720930751902353904327148559471526397")
  const a1206 = BigFloat("-0.861164195027635666673916999665534573351026060987427093314412")
  const a1207 = BigFloat("5.78428813637537220022999785486578436006872789689499172601856")
  const a1208 = BigFloat("3.28807761985103566890460615937314805477268252903342356581925")
  const a1209 = BigFloat("-2.38633905093136384013422325215527866148401465975954104585807")
  const a1210 = BigFloat("-3.25479342483643918654589367587788726747711504674780680269911")
  const a1211 = BigFloat("-2.16343541686422982353954211300054820889678036420109999154887")

  const a1300 = BigFloat("0.895080295771632891049613132336585138148156279241561345991710")
  const a1302 = BigFloat("0.197966831227192369068141770510388793370637287463360401555746")
  const a1303 = BigFloat("-0.0729547847313632629185146671595558023015011608914382961421311")
  const a1305 = BigFloat("-0.851236239662007619739049371445966793289359722875702227166105")
  const a1306 = BigFloat("0.398320112318533301719718614174373643336480918103773904231856")
  const a1307 = BigFloat("3.63937263181035606029412920047090044132027387893977804176229")
  const a1308 = BigFloat("1.54822877039830322365301663075174564919981736348973496313065")
  const a1309 = BigFloat("-2.12221714704053716026062427460427261025318461146260124401561")
  const a1310 = BigFloat("-1.58350398545326172713384349625753212757269188934434237975291")
  const a1311 = BigFloat("-1.71561608285936264922031819751349098912615880827551992973034")
  const a1312 = BigFloat("-0.0244036405750127452135415444412216875465593598370910566069132")

  const a1400 = BigFloat("-0.915176561375291440520015019275342154318951387664369720564660")
  const a1401 = BigFloat("1.45453440217827322805250021715664459117622483736537873607016")
  const a1404 = BigFloat("-0.777333643644968233538931228575302137803351053629547286334469")
  const a1406 = BigFloat("-0.0910895662155176069593203555807484200111889091770101799647985")
  const a1412 = BigFloat("0.0910895662155176069593203555807484200111889091770101799647985")
  const a1413 = BigFloat("0.777333643644968233538931228575302137803351053629547286334469")

  const a1500 = BigFloat(1//10)
  const a1502 = BigFloat("-0.157178665799771163367058998273128921867183754126709419409654")
  const a1514 = BigFloat("0.157178665799771163367058998273128921867183754126709419409654")

  const a1600 = BigFloat("0.181781300700095283888472062582262379650443831463199521664945")
  const a1601 = BigFloat(27//40)
  const a1602 = BigFloat("0.342758159847189839942220553413850871742338734703958919937260")
  const a1604 = BigFloat("0.259111214548322744512977076191767379267783684543182428778156")
  const a1605 = BigFloat("-0.358278966717952089048961276721979397739750634673268802484271")
  const a1606 = BigFloat("-1.04594895940883306095050068756409905131588123172378489286080")
  const a1607 = BigFloat("0.930327845415626983292300564432428777137601651182965794680397")
  const a1608 = BigFloat("1.77950959431708102446142106794824453926275743243327790536000")
  const a1609 = BigFloat(1//10)
  const a1610 = BigFloat("-0.282547569539044081612477785222287276408489375976211189952877")
  const a1611 = BigFloat("-0.159327350119972549169261984373485859278031542127551931461821")
  const a1612 = BigFloat("-0.145515894647001510860991961081084111308650130578626404945571")
  const a1613 = BigFloat("-0.259111214548322744512977076191767379267783684543182428778156")
  const a1614 = BigFloat("-0.342758159847189839942220553413850871742338734703958919937260")
  const a1615 = BigFloat(-27//40)

  const b00 = BigFloat(1//30)
  const b01 = BigFloat(1//40)
  const b02 = BigFloat(1//30)
  const b04 = BigFloat(1//20)
  const b06 = BigFloat(1//25)
  const b08 = BigFloat("0.189237478148923490158306404106012326238162346948625830327194")
  const b09 = BigFloat("0.277429188517743176508360262560654340428504319718040836339472")
  const b10 = BigFloat("0.277429188517743176508360262560654340428504319718040836339472")
  const b11 = BigFloat("0.189237478148923490158306404106012326238162346948625830327194")
  const b12 = BigFloat(-1//25)
  const b13 = BigFloat(-1//20)
  const b14 = BigFloat(-1//30)
  const b15 = BigFloat(-1//40)
  const b16 = BigFloat(1//30)

  const c01 = BigFloat(1//10)
  const c02 = BigFloat("0.539357840802981787532485197881302436857273449701009015505500")
  const c03 = BigFloat("0.809036761204472681298727796821953655285910174551513523258250")
  const c04 = BigFloat("0.309036761204472681298727796821953655285910174551513523258250")
  const c05 = BigFloat("0.981074190219795268254879548310562080489056746118724882027805")
  const c06 = BigFloat(5//6)
  const c07 = BigFloat("0.354017365856802376329264185948796742115824053807373968324184")
  const c08 = BigFloat("0.882527661964732346425501486979669075182867844268052119663791")
  const c09 = BigFloat("0.642615758240322548157075497020439535959501736363212695909875")
  const c10 = BigFloat("0.357384241759677451842924502979560464040498263636787304090125")
  const c11 = BigFloat("0.117472338035267653574498513020330924817132155731947880336209")
  const c12 = BigFloat(5//6)
  const c13 = BigFloat("0.309036761204472681298727796821953655285910174551513523258250")
  const c14 = BigFloat("0.539357840802981787532485197881302436857273449701009015505500")
  const c15 = BigFloat(1//10)
  const c16 = BigFloat(1)

  k00 = h*f(x,         y)
  k01 = h*f(x + c01*h, y + a0100*k00)
  k02 = h*f(x + c02*h, y + a0200*k00 + a0201*k01)
  k03 = h*f(x + c03*h, y + a0300*k00             + a0302*k02)
  k04 = h*f(x + c04*h, y + a0400*k00             + a0402*k02 + a0403*k03)
  k05 = h*f(x + c05*h, y + a0500*k00                         + a0503*k03 + a0504*k04)
  k06 = h*f(x + c06*h, y + a0600*k00                         + a0603*k03 + a0604*k04 + a0605*k05)
  k07 = h*f(x + c07*h, y + a0700*k00                                     + a0704*k04 + a0705*k05 + a0706*k06)
  k08 = h*f(x + c08*h, y + a0800*k00                                                 + a0805*k05 + a0806*k06 + a0807*k07)
  k09 = h*f(x + c09*h, y + a0900*k00                                                 + a0905*k05 + a0906*k06 + a0907*k07 + a0908*k08)
  k10 = h*f(x + c10*h, y + a1000*k00                                                 + a1005*k05 + a1006*k06 + a1007*k07 + a1008*k08 + a1009*k09)
  k11 = h*f(x + c11*h, y + a1100*k00                                                 + a1105*k05 + a1106*k06 + a1107*k07 + a1108*k08 + a1109*k09 + a1110*k10)
  k12 = h*f(x + c12*h, y + a1200*k00                         + a1203*k03 + a1204*k04 + a1205*k05 + a1206*k06 + a1207*k07 + a1208*k08 + a1209*k09 + a1210*k10 + a1211*k11)
  k13 = h*f(x + c13*h, y + a1300*k00             + a1302*k02 + a1303*k03             + a1305*k05 + a1306*k06 + a1307*k07 + a1308*k08 + a1309*k09 + a1310*k10 + a1311*k11 + a1312*k12)
  k14 = h*f(x + c14*h, y + a1400*k00 + a1401*k01                         + a1404*k04             + a1406*k06 +                                                             a1412*k12 + a1413*k13)
  k15 = h*f(x + c15*h, y + a1500*k00             + a1502*k02                                                                                                                                     + a1514*k14)
  k16 = h*f(x + c16*h, y + a1600*k00 + a1601*k01 + a1602*k02             + a1604*k04 + a1605*k05 + a1606*k06 + a1607*k07 + a1608*k08 + a1609*k09 + a1610*k10 + a1611*k11 + a1612*k12 + a1613*k13 + a1614*k14 + a1615*k15)

  # 10th order solution
  sol = y + b00*k00 + b01*k01 + b02*k02 + b04*k04 + b06*k06 + b08*k08 + b09*k09 + b10*k10 + b11*k11 + b12*k12 + b13*k13 + b14*k14 + b15*k15 + b16*k16
  # 8th order error estimate
  err = (k01 - k15) / 360
  return sol, err
end

Answer 1

I have optimized Feagin solvers in DifferentialEquations.jl. You can find the codes starting here. I know your question was about simplifying the code, but after really optimizing them, I believe that is the wrong direction (except for teaching purposes. For teaching, just put A in a sparse matrix, b,c as vectors, and go with it).

The code needs to be almost fully unrolled for a few reasons. First of all, if you test static arrays/matrices for the coefficients of this size (via StaticArrays.jl, they are not good at staying stack allocated. The reason is because they are more larger than what you want to always keep on the stack. This is exacerbated by the fact that A is sparse (even in the lower triangular part that), and so using a full matrix would have large amounts of memory overhead. So you want to keep at least a as stand-alone variables. b and c are easily put into a StaticArray. That simplifies things a little bit. (I tested and it's the same as stand-alone variables. I haven't done this yet on master because it's only v0.5 and I want to keep v0.4 compatibility right now).

But now towards optimization. As @MattB noted, the first thing is to deal with the constants. They should be moved outside of this inner function. Also you should take some pre-caution to make sure that all of the types match (to make the code more general. Maybe you will want to use ArbFloats for better performance?). To add this robustness, you just generate all of the constants inside another function, and you can parse them to a given type T (inside that function, so that way in your main function they are always type T to avoid type instabilities) with

T(parse(BigFloat,"0.105352113571753019691496032887878162227673083080523884041670"))

Just moving the constants out of the loop should gives you orders of magnitude speedup.

If your y is 1-dimension (i.e. just a number), then the last optimization I could find is to group the variables in long summations by 4's. For example, (a+b+c+d)+(e+f+g+h)+... instead of a+b+c+d+e+f+g+h+... The reason is due to a compiler inefficiency. Putting these parenthesis around can increase the speed of the code by a few times.

A bit more fine-tuned of an optimization is to now reduce the temporaries. Anytime you have a vectorized operation like A+B it stores the result in a temporary variable. So in the case where y and all of the k's are vectors, there are a lot of temporaries that are allocated each step. First you can get rid of some by making f an in-place function which takes in a (du,u,t) and doesn't return the value, but instead writes it into du. In the spot of the du you will put the k's. Once that's in-place, then you will want to de-vectorize all of the large summations. I created a Gist where I tested which is the fastest way to do this. The winner is to do the long-single loop with parenthesis grouping the terms by 4. You'd want to do that to every line of a's in the k = f lines, and you'd want to do that in the solution update part. You can re-use the same du in every line because it's only used for one part, and thus you only need a few cache variables around for the whole approximation.

The result of that optimization is an utter mess. Refer to the DifferentialEquations.jl source code again. However timings show that even for small y vectors you can get about 1/3 faster. Then you can change those cache arrays to StaticArrays to probably get major speed gains (but require v0.5).

So the resulting code is no simpler, but it's much faster! Note that DifferentialEquations.jl has these hand-unrolled devecotorized implementations of also the order 12 and 14 methods if you'd like to take a look. It gets quite tedious to do, but sometimes you have to do it if you want the most optimal code.