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
let
scope (the tricky part is that you need to declare your function as "global" to escape the scope). See, for example, hownpartitions
in base caches its results. \$\endgroup\$ – Matt B. Apr 16 '15 at 15:14