cowell4.lua

cowell4.lua

-- Lua による常微分方程式の数値解法 / 4 段のカウエル法
function cowell4(funcs, t0, inits, h)
   local unpack = unpack or table.unpack
  
   local t0    = t0
   local inits = inits

   local function maxOfErr(listA, listB)
            local sute = {}
            for i = 1, #listA do
               sute[i] = math.abs(listA[i] - listB[i])
            end
            return math.max(unpack(sute))
         end
   local function staticVals(funcs, inits, h)
            local dim  = #funcs
            local w    = {}
            local f    = {{},{},{},{}}
            local temp = {}
            for i = 1, 4 do
               for j = 1, dim do
                  f[i][j] = funcs[j](unpack(inits[i]))
               end
            end  
            for j = 1, dim do
               w[j]    = 2 * inits[4][j] - inits[3][j] + h*h * (204 * f[4][j] + 14 * f[3][j] + 4 * f[2][j] - f[1][j]) / 240
               temp[j] = 2 * f[4][j] - f[3][j]
            end
            return w, temp -- {変化しない部分}, {假定値}
         end
   local function approx(funcs, w, temp, h) -- ({微分方程式}, {変化しない部分}, {假定値}, 刻み幅)
            local dim     = #funcs
            local approx1 = {}
            local approx2 = {}
            local f       = {}
            for j = 1, dim do
               approx1[j] = w[j] + h*h * 19 * temp[j]/240 -- 第 1 近似値
            end
            for j = 1, dim do
               f[j] = funcs[j](unpack(approx1)) 
            end
            for j = 1, dim do
               approx2[j] = w[j] + h*h * 19 * f[j]/240 -- 第 2 近似値
            end
            return approx1, approx2, maxOfErr(approx1, approx2) -- {第 1 近似値}, {第 2 近似値}, 最大誤差  
         end
   local function approx2(funcs, w, preApprox, h) -- ({微分方程式}, {変化しない部分}, {前の近似値}, 刻み幅)
            local dim    = #funcs
            local f      = {}
            local approx = {}
            for j = 1, dim do
               f[j] = funcs[j](unpack(preApprox)) 
            end
            for j = 1, dim do
               approx[j] = w[j] + h*h * 19 * f[j]/240 -- 次の近似値
            end
            return approx, maxOfErr(preApprox, approx) -- {次の近似値}, (前の近似値と次の近似値との最大誤差)  
         end

   return function()
              local count = 0 
              local w, temp = staticVals(funcs, inits, h)                 -- まず変化しない部分と最初の假定値とを求めて、
              local approx1st, approx2nd, err = approx(funcs, w, temp, h) -- 最初の第 1 近似値、第 2 近似値、両者の差を求めて、
              while(err > 1e-15) do                                       -- 差が存在しなかったらループを飛ばして初期値を更新するが、差が存在していたらループに入って、
                  approx2nd, err = approx2(funcs, w, approx2nd, h)        -- 差がなくなるまで近似値の計算を繰り返して、差がなくなったらループを抜けて初期値を更新する。
                  count = count + 1
              end
              table.remove(inits, 1)
              inits[4] = approx2nd
              t0 = t0 + h
              return t0, inits[1], count -- 独立変数の値, {4 組の初期値の先頭の初期値}, 近似回数  という形で返す。
           end
end


-- 確かめる
do

-- この聯立二階微分方程式で確かめる。
local function xDotDot(x, y) return -(2.95912208e-4) * x / (math.sqrt(x^2 + y^2)^3) end
local function yDotDot(x, y) return -(2.95912208e-4) * y / (math.sqrt(x^2 + y^2)^3) end

-- 4 組の初期値にはこの値を使う。 {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}
local inits = {{1.0989720 , 0        },
               {1.0509145 , 0.4038387},
               {0.9168459 , 0.7754143},
               {0.7200885 , 1.0952558}}
local h  = 20
local t0 = 0 

-- 1 ステップだけ計算する函数を作って、
local a = cowell4({xDotDot, yDotDot}, t0, inits, h)

-- 必要な回数だけその函数を実行する。
for i = 1, 80 do 
   local t, inits, count = a()
   print(t, table.concat(inits, ", "), count) 
end

end