cowell7.lua

cowell7.lua

-- Lua で常微分方程式を解く / 7 段のカウェル法
function cowell7(funcs, t0, inits, h)
   local unpack = unpack or table.unpack
   local t0     = t0
   local inits  = inits

   function staticVals(funcs, inits, h)
      local dim  = #funcs
      local w    = {}
      local f    = {{},{},{},{},{},{},{}}
      local temp = {}
      for i = 1, 7 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[7][j] - inits[6][j] + h*h * (55324 * f[7][j] - 6297 * f[6][j] + 14598 * f[5][j] - 11477 * f[4][j] + 5568 * f[3][j] - 1551 * f[2][j] + 190 * f[1][j]) / 60480
         temp[j] = 2 * f[7][j] - f[6][j]
      end
      return w, temp -- {変化しない部分}, {假定値}
   end
   function approx(funcs, w, temp, h) -- ({微分方程式}, {変化しない部分}, {假定値}, 刻み幅)
      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 dim     = #funcs
      local approx1 = {}
      local approx2 = {}
      local f       = {}
      for j = 1, dim do
         approx1[j] = w[j] + h*h * 4125 * temp[j]/60480 -- 第 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 * 4125 * f[j]/60480 -- 第 2 近似値
      end
      return approx1, approx2, maxOfErr(approx1, approx2) -- {第 1 近似値}, {第 2 近似値}, 最大誤差  
   end
   function APPROX(funcs, w, preApprox, h) -- ({微分方程式}, {変化しない部分}, {前の近似値}, 刻み幅)
      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 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 * 4125 * f[j]/60480 -- 次の近似値
      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 = APPROX(funcs, w, approx2nd, h)         -- 差がなくなるまで近似値の計算を繰り返して、差がなくなったらループを抜けて初期値を更新する。
                  count = count + 1
              end
              table.remove(inits, 1)
              inits[7] = approx2nd
              t0 = t0 + h
              return t0, inits[1], count -- 独立変数の値, {7 組の初期値の先頭の初期値}, 近似回数
           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

-- 7 組の初期値にはこの値を使う。 {{x1, y1}, {x2, y2}, ..., {x6, y6}, {x7, y7}}
local inits = {{ 1.0989720  , 0.0      },
               { 1.0509145  , 0.4038387},
               { 0.9168459  , 0.7754143},
               { 0.7200885  , 1.0952558},
               { 0.4849475  , 1.3582345},
               { 0.2301295  , 1.5678560},
               {-0.03197416 , 1.7309904}}
local h  = 20
local t0 = 0

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

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


end