TI-Nspire & Lua / Shanks のルンゲクッタ法公式

参考: 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

f:id:ti-nspire:20170508081417p:plain:w400
f:id:ti-nspire:20170508081456p:plain:w400