参考: Shanks, E.B., (1965), Higher Order Approximations of Runge-Kutta Type, NASA TN D-2920 (https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19650022581.pdf)
function rkShanks(funcs, t0, inits, h, numOfDiv) local unpack = unpack or table.unpack local dim = #funcs local numOfDiv = numOfDiv or 1 local h = h / numOfDiv local fnc01, fnc02, fnc03, fnc04, fnc05, fnc06, fnc07, fnc08, fnc09, fnc10, fnc11, fnc12 = {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {} local tmp01, tmp02, tmp03, tmp04, tmp05, tmp06, tmp07, tmp08, tmp09, tmp10, tmp11 = {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {} local t0 = t0 local inits = inits return function(order) local order = order or 812 -- order は 44、55、66、77、79、810、812 のいずれかで指定する。デフォルトは 812。 for i = 1, numOfDiv do if order == 44 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h/100 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h/100, unpack(tmp01)) ; tmp02[j] = inits[j] + h/245 * (-4278 * fnc01[j] + 4425 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h*3/5, unpack(tmp02)) ; tmp03[j] = inits[j] + h/8791 * (524746 * fnc01[j] - 532125 * fnc02[j] + 16170 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h , unpack(tmp03)) end for j = 1, dim do inits[j] = inits[j] + h/70092 * (-179124 * fnc01[j] + 200000 * fnc02[j] + 40425 * fnc03[j] + 8791 * fnc04[j]) end elseif order == 55 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h/9000 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h/9000, unpack(tmp01)) ; tmp02[j] = inits[j] + h/10 * (-4047 * fnc01[j] + 4050 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h*3/10, unpack(tmp02)) ; tmp03[j] = inits[j] + h/8 * (20241 * fnc01[j] - 20250 * fnc02[j] + 15 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h*3/4 , unpack(tmp03)) ; tmp04[j] = inits[j] + h/81 * (-931041 * fnc01[j] + 931500 * fnc02[j] - 490 * fnc03[j] + 112 * fnc04[j]) end for j = 1, dim do fnc05[j] = funcs[j](t0 + h , unpack(tmp04)) end for j = 1, dim do inits[j] = inits[j] + h/1134 * (105 * fnc01[j] + 500 * fnc03[j] + 448 * fnc04[j] + 81 * fnc05[j]) end elseif order == 66 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h/300 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h/300 , unpack(tmp01)) ; tmp02[j] = inits[j] + h/5 * (-29 * fnc01[j] + 30 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h/5 , unpack(tmp02)) ; tmp03[j] = inits[j] + h/5 * (323 * fnc01[j] - 330 * fnc02[j] + 10 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h*3/5 , unpack(tmp03)) ; tmp04[j] = inits[j] + h/810 * (-510104 * fnc01[j] + 521640 * fnc02[j] - 12705 * fnc03[j] + 1925 * fnc04[j]) end for j = 1, dim do fnc05[j] = funcs[j](t0 + h*14/15, unpack(tmp04)) ; tmp05[j] = inits[j] + h/77 * (-417923 * fnc01[j] + 427350 * fnc02[j] - 10605 * fnc03[j] + 1309 * fnc04[j] - 54 * fnc05[j]) end for j = 1, dim do fnc06[j] = funcs[j](t0 + h , unpack(tmp05)) end for j = 1, dim do inits[j] = inits[j] + h/3696 * (198 * fnc01[j] + 1225 * fnc03[j] + 1540 * fnc04[j] + 810 * fnc05[j] -77 * fnc06[j]) end elseif order == 77 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h/192 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h/192, unpack(tmp01)) ; tmp02[j] = inits[j] + h/6 * (-15 * fnc01[j] + 16 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h/6 , unpack(tmp02)) ; tmp03[j] = inits[j] + h/186 * (4867 * fnc01[j] - 5072 * fnc02[j] + 298 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h/2 , unpack(tmp03)) ; tmp04[j] = inits[j] + h/31 * (-19995 * fnc01[j] + 20896 * fnc02[j] - 1025 * fnc03[j] + 155 * fnc04[j]) end for j = 1, dim do fnc05[j] = funcs[j](t0 + h , unpack(tmp04)) ; tmp05[j] = inits[j] + h/5022 * (-469805 * fnc01[j] + 490960 * fnc02[j] - 22736 * fnc03[j] + 5580 * fnc04[j] + 186 * fnc05[j]) end for j = 1, dim do fnc06[j] = funcs[j](t0 + h*5/6, unpack(tmp05)) ; tmp06[j] = inits[j] + h/2604 * (914314 * fnc01[j] - 955136 * fnc02[j] + 47983 * fnc03[j] - 6510 * fnc04[j] - 558 * fnc05[j] + 2511 * fnc06[j]) end for j = 1, dim do fnc07[j] = funcs[j](t0 + h , unpack(tmp06)) end for j = 1, dim do inits[j] = inits[j] + h/300 * (14 * fnc01[j] + 81 * fnc03[j] + 110 * fnc04[j] + 81 * fnc06[j] + 14 * fnc07[j]) end elseif order == 79 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h*2/9 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h*2/9, unpack(tmp01)) ; tmp02[j] = inits[j] + h/12 * ( fnc01[j] + 3 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h/3 , unpack(tmp02)) ; tmp03[j] = inits[j] + h/8 * ( fnc01[j] + 3 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h/2 , unpack(tmp03)) ; tmp04[j] = inits[j] + h/216 * (23 * fnc01[j] + 21 * fnc03[j] - 8 * fnc04[j]) end for j = 1, dim do fnc05[j] = funcs[j](t0 + h/6 , unpack(tmp04)) ; tmp05[j] = inits[j] + h/729 * (-4136 * fnc01[j] - 13584 * fnc03[j] + 5264 * fnc04[j] + 13104 * fnc05[j]) end for j = 1, dim do fnc06[j] = funcs[j](t0 + h*8/9, unpack(tmp05)) ; tmp06[j] = inits[j] + h/151632 * (105131 * fnc01[j] + 302016 * fnc03[j] - 107744 * fnc04[j] - 284256 * fnc05[j] + 1701 * fnc06[j]) end for j = 1, dim do fnc07[j] = funcs[j](t0 + h/9 , unpack(tmp06)) ; tmp07[j] = inits[j] + h/1375920 * (-775229 * fnc01[j] - 2770950 * fnc03[j] + 1735136 * fnc04[j] + 2547216 * fnc05[j] + 81891 * fnc06[j] + 328536 * fnc07[j]) end for j = 1, dim do fnc08[j] = funcs[j](t0 + h*5/6, unpack(tmp07)) ; tmp08[j] = inits[j] + h/251888 * (23569 * fnc01[j] - 122304 * fnc03[j] - 20384 * fnc04[j] + 695520 * fnc05[j] - 99873 * fnc06[j] - 466560 * fnc07[j] + 241920 * fnc08[j]) end for j = 1, dim do fnc09[j] = funcs[j](t0 + h , unpack(tmp08)) end for j = 1, dim do inits[j] = inits[j] + h/2140320 * (110201 * fnc01[j] + 767936 * fnc04[j] + 635040 * fnc05[j] - 59049 * fnc06[j] - 59049 * fnc07[j] + 635040 * fnc08[j] + 110201 * fnc09[j]) end elseif order == 810 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h*4/27 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h*4/27, unpack(tmp01)) ; tmp02[j] = inits[j] + h/18 * ( fnc01[j] + 3 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h*2/9 , unpack(tmp02)) ; tmp03[j] = inits[j] + h/12 * ( fnc01[j] + 3 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h/3 , unpack(tmp03)) ; tmp04[j] = inits[j] + h/8 * ( fnc01[j] + 3 * fnc04[j]) end for j = 1, dim do fnc05[j] = funcs[j](t0 + h/2 , unpack(tmp04)) ; tmp05[j] = inits[j] + h/54 * (13 * fnc01[j] - 27 * fnc03[j] + 42 * fnc04[j] + 8 * fnc05[j]) end for j = 1, dim do fnc06[j] = funcs[j](t0 + h*2/3 , unpack(tmp05)) ; tmp06[j] = inits[j] + h/4320 * (389 * fnc01[j] - 54 * fnc03[j] + 966 * fnc04[j] - 824 * fnc05[j] + 243 * fnc06[j]) end for j = 1, dim do fnc07[j] = funcs[j](t0 + h/6 , unpack(tmp06)) ; tmp07[j] = inits[j] + h/20 * (-231 * fnc01[j] + 81 * fnc03[j] - 1164 * fnc04[j] + 656 * fnc05[j] - 122 * fnc06[j] + 800 * fnc07[j]) end for j = 1, dim do fnc08[j] = funcs[j](t0 + h , unpack(tmp07)) ; tmp08[j] = inits[j] + h/288 * (-127 * fnc01[j] + 18 * fnc03[j] - 678 * fnc04[j] + 456 * fnc05[j] - 9 * fnc06[j] + 576 * fnc07[j] + 4 * fnc08[j]) end for j = 1, dim do fnc09[j] = funcs[j](t0 + h*5/6 , unpack(tmp08)) ; tmp09[j] = inits[j] + h/820 * (1481 * fnc01[j] - 81 * fnc03[j] + 7104 * fnc04[j] - 3376 * fnc05[j] + 72 * fnc06[j] - 5040 * fnc07[j] - 60 * fnc08[j] + 720 * fnc09[j]) end for j = 1, dim do fnc10[j] = funcs[j](t0 + h , unpack(tmp09)) end for j = 1, dim do inits[j] = inits[j] + h/840 * (41 * fnc01[j] + 27 * fnc04[j] + 272 * fnc05[j] + 27 * fnc06[j] + 216 * fnc07[j] + 216 * fnc09[j] + 41 * fnc10[j]) end elseif order == 812 then for j = 1, dim do fnc01[j] = funcs[j](t0 , unpack(inits)) ; tmp01[j] = inits[j] + h/9 * ( fnc01[j]) end for j = 1, dim do fnc02[j] = funcs[j](t0 + h/9 , unpack(tmp01)) ; tmp02[j] = inits[j] + h/24 * ( fnc01[j] + 3 * fnc02[j]) end for j = 1, dim do fnc03[j] = funcs[j](t0 + h/6 , unpack(tmp02)) ; tmp03[j] = inits[j] + h/16 * ( fnc01[j] + 3 * fnc03[j]) end for j = 1, dim do fnc04[j] = funcs[j](t0 + h/4 , unpack(tmp03)) ; tmp04[j] = inits[j] + h/500 * (29 * fnc01[j] + 33 * fnc03[j] - 12 * fnc04[j]) end for j = 1, dim do fnc05[j] = funcs[j](t0 + h/10 , unpack(tmp04)) ; tmp05[j] = inits[j] + h/972 * (33 * fnc01[j] + 4 * fnc04[j] + 125 * fnc05[j]) end for j = 1, dim do fnc06[j] = funcs[j](t0 + h/6 , unpack(tmp05)) ; tmp06[j] = inits[j] + h/36 * (-21 * fnc01[j] + 76 * fnc04[j] + 125 * fnc05[j] - 162 * fnc06[j]) end for j = 1, dim do fnc07[j] = funcs[j](t0 + h/2 , unpack(tmp06)) ; tmp07[j] = inits[j] + h/243 * (-30 * fnc01[j] - 32 * fnc04[j] + 125 * fnc05[j] + 99 * fnc07[j]) end for j = 1, dim do fnc08[j] = funcs[j](t0 + h*2/3, unpack(tmp07)) ; tmp08[j] = inits[j] + h/324 * (1175 * fnc01[j] - 3456 * fnc04[j] - 6250 * fnc05[j] + 8424 * fnc06[j] + 242 * fnc07[j] - 27 * fnc08[j]) end for j = 1, dim do fnc09[j] = funcs[j](t0 + h/3 , unpack(tmp08)) ; tmp09[j] = inits[j] + h/324 * (293 * fnc01[j] - 852 * fnc04[j] - 1375 * fnc05[j] + 1836 * fnc06[j] - 118 * fnc07[j] + 162 * fnc08[j] + 324 * fnc09[j]) end for j = 1, dim do fnc10[j] = funcs[j](t0 + h*5/6, unpack(tmp09)) ; tmp10[j] = inits[j] + h/1620 * (1303 * fnc01[j] - 4260 * fnc04[j] - 6875 * fnc05[j] + 9990 * fnc06[j] + 1030 * fnc07[j] + 162 * fnc10[j]) end for j = 1, dim do fnc11[j] = funcs[j](t0 + h*5/6, unpack(tmp10)) ; tmp11[j] = inits[j] + h/4428 * (-8595 * fnc01[j] + 30720 * fnc04[j] + 48750 * fnc05[j] - 66096 * fnc06[j] + 378 * fnc07[j] - 729 * fnc08[j] - 1944 * fnc09[j] - 1296 * fnc10[j] + 3240 * fnc11[j]) end for j = 1, dim do fnc12[j] = funcs[j](t0 + h , unpack(tmp11)) end for j = 1, dim do inits[j] = inits[j] + h/840 * (41 * fnc01[j] + 216 * fnc06[j] + 272 * fnc07[j] + 27 * fnc08[j] + 27 * fnc09[j] + 36 * fnc10[j] + 180 * fnc11[j] + 41 * fnc12[j]) end end t0 = t0 + h end return t0, unpack(inits) end end --- 確かめる。 do -- この聯立微分方程式を解く。 local function xDot(t, x, y) return y end local function yDot(t, x, y) return t - x end -- 1 ステップだけ計算する函数を作って、 local ShanksA = rkShanks({xDot, yDot}, 0, {0, 0}, 0.125, _) local ShanksB = rkShanks({xDot, yDot}, 0, {0, 0}, 0.125, _) -- その函数を必要な回数だけ実行する。 for i = 1, 7.5/0.125 do print(ShanksA(44)) -- 4 次 4 段と print(ShanksB(812)) -- 8 次 12 段とを比べる。 end end