TI-Nspire & Lua / Fehlberg 法 6 of 6 / ローレンツアトラクターを描く / 1 個の解の軌跡を描く

xz 面の解の軌跡を描く。y の値は圓の大きさで表現する。

function fehlberg(funcs, t0, inits, h, tol)
   local unpack = unpack or table.unpack
   local Ln     = math.log
   local Floor  = math.floor
   local Abs    = math.abs
   local Max    = math.max  
   local Pow    = math.pow

   local dim = #funcs

   local t0    = t0
   local inits = inits

   local function floorB(num)
            if (num > 0) and (num < 1) then 
               return Pow(2, Floor(Ln(num)/0.69314718055995))
            else 
               return Floor(num)
            end
         end
   local function hNew(h, err, tol)
            if err > tol then
               return 0.9 * h * Pow(tol/err, 0.25)
            else
               return h
            end
         end
   local function maxOfErr(listA, listB)
            local sute = {}
            for i = 1, #listA do
               sute[i] = Abs(listA[i] - listB[i])
            end
            return Max(unpack(sute))
         end
   local function preCalc(numOfDiv)
            local f1  , f2  , f3  , f4  , f5  , f6 = {}, {}, {}, {}, {}, {}
            local tmp1, tmp2, tmp3, tmp4, tmp5     = {}, {}, {}, {}, {}

            local inits4   = {}
            local preInits = {unpack(inits)}
            local preT0    = t0
            local preH     = h/numOfDiv

            for i = 1, numOfDiv do
               for j = 1, dim do f1[j] = funcs[j](preT0                            , unpack(preInits)); tmp1[j] = preInits[j] + preH * (0.25)               * f1[j]                                                                                                                  end 
               for j = 1, dim do f2[j] = funcs[j](preT0 + preH * (0.25)            , unpack(tmp1))    ; tmp2[j] = preInits[j] + preH * (0.03125)            * (3                  * f1[j] + 9    * f2[j])                                                                            end
               for j = 1, dim do f3[j] = funcs[j](preT0 + preH * (0.375)           , unpack(tmp2))    ; tmp3[j] = preInits[j] + preH * (4.5516613563951e-4) * (1932               * f1[j] - 7200 * f2[j] + 7296              * f3[j])                                                end
               for j = 1, dim do f4[j] = funcs[j](preT0 + preH * (0.92307692307692), unpack(tmp3))    ; tmp4[j] = preInits[j] + preH                        * ((2.0324074074074)  * f1[j] - 8    * f2[j] + (7.1734892787524) * f3[j] - (0.20589668615984) * f4[j])                   end
               for j = 1, dim do f5[j] = funcs[j](preT0 + preH                     , unpack(tmp4))    ; tmp5[j] = preInits[j] + preH                        * ((-0.2962962962963) * f1[j] + 2    * f2[j] - (1.3816764132554) * f3[j] + (0.45297270955166) * f4[j] - (0.275) * f5[j]) end
               for j = 1, dim do f6[j] = funcs[j](preT0 + preH * (0.5)             , unpack(tmp5))                                                                                                                                                                                   end

               if numOfDiv == 1 then -- 内部分割数が 1 のときだけ 4 次の解を計算する。
                  for j = 1, dim do   inits4[j] = preInits[j] + preH * ((0.11574074074074) * f1[j] + (0.54892787524366)  * f3[j] + (0.53533138401559)   * f4[j] - (0.2)  * f5[j]) end
               end
               for j = 1, dim do preInits[j] = preInits[j] + preH * ((0.11851851851852) * f1[j] + (0.51898635477583) * f3[j] + (0.50613149034202) * f4[j] - (0.18) * f5[j] + (0.036363636363636) * f6[j]) end

               preT0 = preT0 + preH
            end
      
      return preT0, inits4, preInits 
   end
   
   return function()
             local numOfDiv = 1
             local t, a4, a5 = preCalc(numOfDiv)       -- とりあえず元の刻み幅で 1 回だけ計算し、
             local err = maxOfErr(a4, a5)  
             if err < tol then                         -- 誤差が許容範囲内であったら、
                t0, inits = t, a5                      -- 初期値を更新するが、
             else                                      -- 誤差が許容範囲外であったら、
                numOfDiv = h/floorB(hNew(h, err, tol)) -- 新しい内部分割数を求めて、
                t, a4, a5  = preCalc(numOfDiv)         -- その内部分割数で計算し直し、
                t0, inits  = t, a5                     -- 初期値を更新する。
             end
             return t0, inits, numOfDiv
          end
   
end


---------------------------------------------------------------------------------
-- ローレンツアトラクターで確かめる。
do
require "color"
platform.window:setBackgroundColor(color.black)

local W, H   = platform.window:width(), platform.window:height()
local XC, YC = W/2, H
local UNIT   = W/73

local Max         = math.max
local tableRemove = table.remove
 
local function fillCircle(centerX, centerY, radius, color, gc)
         local diameter = Max(0, radius) * 2
         gc:setColorRGB(color)
         gc:fillArc(centerX - radius, centerY - radius, diameter, diameter, 0, 360)
      end
local function xy2screen(x, y, xc, yc, unit)
         return xc + x * unit, yc - y * unit  
      end

-- この微分方程式を解く。
local function xDot(t, x, y, z) return 10 * (y - x)                   end
local function yDot(t, x, y, z) return  x * (28 - z) - y              end
local function zDot(t, x, y, z) return  x * y - (2.6666666666667) * z end

-- この初期値で解く。
local t0         = 0
local x0, y0, z0 = 1, 0, 0
local h          = 1/(2^6)
local xz         = {xy2screen(x0, z0, XC, YC, UNIT)}
local t, div = t0, h/h
local a5 = {}

-- 1 ステップだけ計算する函数を作って、
local fe = fehlberg({xDot, yDot, zDot}, t0, {x0, y0, z0}, h, 1e-10)

function on.timer()
   t, a5, div = fe()
   xz[#xz+1], xz[#xz+2] = xy2screen(a5[1], a5[3], XC, YC, UNIT)
   if #xz > 1000 then
      tableRemove(xz, 1) ; tableRemove(xz, 1)
   end
   --print(t, "{"..table.concat(a5, ", ").."}", div, table.concat(xz, " , "))
   platform.window:invalidate()
end
function on.enterKey()
   timer.start(0.01)
end
function on.paint(gc) 
   fillCircle(xz[#xz-1], xz[#xz], ((a5[2] or 0) + 10) * 1.5, 0x00E5EE, gc)
   gc:setColorRGB(0x00BFFF)
   gc:drawString("t : "       ..(t or "--"  ), 10, 0)
   gc:drawString("numOfDiv : "..(div or "--"), 10, 20)
   gc:drawPolyLine(xz)
end

end