常微分方程式の数値解法

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 円制限三体問題 2 / Graphs アプリケーションで動かしてみる / Fehlberg 法

ここでも、L4 を原点に置き、太陽を廻轉中心とする廻轉座標に描く。 Fehlberg = class() function Fehlberg:init(funcs, t0, inits, h, tol) self.funcs = funcs self.t0 = t0 self.inits = inits self.h = h self.step = h self.numOfDiv = nil self.dim = …

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 円制限三体問題 1 / トロヤ群小惑星の秤動 / Fehlberg 法

参考: パソコンで見る天体の動き, p.162-166 Excelで操る! ここまでできる科学技術計算, p.80-86 L4 を原点に置き、太陽と木星とを結ぶ直線を x 軸に平行に固定して (つまり廻轉座標系に) 描く。ここでは両テキストと同じく Fehlberg 法を使ってみる。L4 か…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / ばね振り子 2 of 2 / 軌跡を描く / 補外法

Extrapo = class() function Extrapo:init(funcs, t0, inits, h, numOfDiv) self.Sin = math.sin self.Cos = math.cos self.historicData = {} local unpack = unpack or table.unpack self.numOfDiv = numOfDiv or 1 self.funcs = funcs self.step = h self…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / ばね振り子 1 / 補外法

参考: ばね振り子 - epii's physics notes Extrapo = class() function Extrapo:init(funcs, t0, inits, h, numOfDiv) self.Sin = math.sin self.Cos = math.cos local unpack = unpack or table.unpack self.numOfDiv = numOfDiv or 1 self.funcs = funcs s…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 5 of 5 / 補外法 / 質点同士の接近時のみ正則化する / それに意味のないとき

参考: パソコンで見る天体の動き, p.166-175原点に固定してある超重量物との距離が 0.65 よりも近くなったときだけ假想時間軸に移る。 TwoBodies = class() function TwoBodies:init(funcs, t0, inits, h, numOfDiv, inWhichDomain) self.unpack = unpack or…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 4 / 補外法 / 正則化した方程式を解く / それに意味のないとき

参考: パソコンで見る天体の動き, p.166-175 Levi-Civita の変換で正則化した方程式を解く。假想時間軸の進みは一様であるが、現実時間軸の進みは、質点同士 (一方の質点は原点に固定してある) の近づいたときに遅くなる。誤差を見積もらない限り微分方程式…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 3 / Fehlberg 法 / 刻み幅を細かく、精度を厳しくする / それに意味のないとき

時間を細かく刻めば改善はする。ここでは Fehlberg 法を使っているので、精度を厳しくしても改善する。しかしいずれも彌縫策に過ぎない。 TwoBodies = class() function TwoBodies:init(funcs, t0, inits, h, tol) self.step = h self.funcs = funcs self.t0…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 2 / Fehlberg 法 / それに意味のないとき / 質点同士が極端に近づくとき

参考: パソコンで見る天体の動き, p.166 簡単にするため、質点同士の質量差が極端であるものとし、極端に重いほうの質点 (たとえば太陽とする) が原点から動かないものとする。長半径 1、離心率 0.999、周期 2 の条件で極端に軽いほうの質点 (何らかの天体と…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 体問題 1 / Fehlberg 法 / それに意味のないとき

Lua の実行画面は使わずに Graphs アプリケーションに描画してみる。再生、一時停止、リセットの各操作はスライダーをボタン代わりにして行う。 TwoBodies = class() function TwoBodies:init(funcs, t0, inits, h, tol) self.step = h self.funcs = funcs s…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 聯成振動 2 / Nyström 法 / ゴムひもの伸び縮みを表現する

CoupledVib = class() function CoupledVib:init(funcs, t0, inits, initsDot, h, numOfDiv) self.fillLine = (function() local POW = math.pow local ATAN = math.atan2 local MIN = math.min local SQRT = math.sqrt local SIN = math.sin local COS = ma…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 聯成振動 1 / Nyström 法 / それに意味のないとき

参考: Excelで学ぶ微分方程式, p.204-209テキストは古典的 Runge-Kutta 法を使っているが、ここでは、特殊な方程式に対する Nyström 法を使ってみる。 縦の点線が平衡点: CoupledVib = class() function CoupledVib:init(funcs, t0, inits, initsDot, h, num…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 重振子でバタフライ効果を見る / 古典的 Runge-Kutta 法

L2 の長さを 0.01% だけ違えた 2 組の 2 重振子を同一座標に描く。t = 6 あたりまでは同じ軌道を描くが、t = 10 を過ぎると全然違う軌道になる。 DoublePend = class() function DoublePend:init(funcs, t0, inits, h, numOfDiv, len1, radius1, len2, radiu…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 2 重振子 / 古典的 Runge-Kutta 法

参考: 工業力学 (機械工学基礎講座), p.229-231 DoublePend = class() function DoublePend:init(funcs, t0, inits, h, numOfDiv, len1, radius1, len2, radius2) self.len1 = len1 self.radius1 = radius1 self.len2 = len2 self.radius2 = radius2 self.hi…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / Pendulum Waves / 古典的 Runge-Kutta 法

長さの違う複数の振子を同じ初期角から一斉に振る。各糸の長さは、たとえば同じ 30 秒間に振子 1 は 15 往復、振子 2 は 16 往復、振子 3 は 17 往復 ...... するように設定する。参考: Pendulum Waves Pend = class() function Pend:init(funcs, t0, inits,…

TI-Nspire & Lua / 常微分方程式の数値解法 / 応用例 / 単振子 / 古典的 Runge-Kutta 法

参考: Excelで学ぶ微分方程式, p.171-176 Pend = class() function Pend:init(funcs, t0, inits, h, numOfDiv, len, radius) self.len = len self.radius = radius self.Unpack = unpack or table.unpack self.Sin = math.sin self.Cos = math.cos self.func…

TI-Nspire & Lua / 常微分方程式の数値解法のまとめ (『パソコンで見る天体の動き』) / クラス化する

同じ処理を今度はクラス化する。 1. 古典的 Runge-Kutta 法 構文: rkClassic({funcs}, t0, {inits}, h [, numOfDiv])▶ソースコードを表示する/非表示にする rkClassicClass.lua 実行結果 (左から t, x, y) 2. Shanks による 12 段 8 次の Runge-Kutta 法 構…

TI-Nspire & Lua / 常微分方程式の数値解法のまとめ (『パソコンで見る天体の動き』) / クロージャを利用する

『パソコンで見る天体の動き』に解説してある常微分方程式の数値解法を Euler 法以外全部 Lua で実装したのでここにまとめておく。全部クロージャを利用している。 ここで言う「特殊な方程式」とは独立変数も 1 次微分の項も含まない 2 階微分方程式のことで…

TI-Nspire & Lua / 特殊な方程式に対する補外法

参考: パソコンで見る天体の動き, p.131-142 function gragg(funcs, t0, inits, initsDot, h) local unpack = unpack or table.unpack local t0 = t0 local inits = inits local initsDot = initsDot local function list2oneColMat(list) local oneColMat =…

TI-Nspire & Lua / 補外法 6 of 6 / まとめ

参考: パソコンで見る天体の動き, p.131-142 function extrapolation(funcs, t0, inits, h) local unpack = unpack or table.unpack local t0 = t0 local inits = inits local function list2oneColMat(list) local oneColMat = {} for i = 1, #list do oneC…

TI-Nspire & Lua / 補外法 5 / 補外にかける点列を中点法で求める

参考: パソコンで見る天体の動き, p.131-142 function newMat(numRows, numCols, val) local mat = {} for i = 1, numRows do mat[i] = {} for j = 1, numCols do mat[i][j] = val or nil end end return mat end function midPoints(funcs, t0, inits, h, s…

TI-Nspire & Lua / 補外法 4 / 何かを偶数次多項式で補外する / Neville の算法で補外値だけを求める / 刻み幅を順次 1/2 にした場合 / 複数の点列に適用する

参考: パソコンで見る天体の動き, p.131-142 local function list2oneColMat(list) local oneColMat = {} for i = 1, #list do oneColMat[i] = {} oneColMat[i][1] = list[i] end return oneColMat end local function neville(...) local listOfLists = {..…

TI-Nspire & Lua / 補外法 3 / 何かを偶数次多項式で補外する / Neville の算法で補外値だけを求める / 刻み幅を順次 1/2 にした場合

参考: パソコンで見る天体の動き, p.131-142 h = {1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128} a = {{0.53627861}, {0.53951755}, {0.54011910}, {0.54025730}, {0.54029111}, {0.54029951}, {0.54030161}} for i = 2, #h do for j = 2, i do --a[i][j] = a[i]…

TI-Nspire & Lua / 補外法 2 / 何かを偶数次多項式で補外する / Neville の算法で補外値だけを求める

参考: パソコンで見る天体の動き, p.131-142既知のポイントの数列は前回と同じとする。 h = {4, 3, 2, 1, 1/2} a = {{277}, {-472}, {-167}, {-8}, {277/64}} for i = 2, #h do for j = 2, i do a[i][j] = a[i][j-1] + (a[i][j-1] - a[i-1][j-1]) / ((h[i-j+…

TI-Nspire & Lua / 補外法 1 / 何かを偶数次多項式で補外する / 聯立方程式を解いて求める

参考: パソコンで見る天体の動き, p.131-142既知のいくつかのポイントを通過する未知の曲線が偶数次多項式であると假定してその既知のポイント群の外側に存在する未知のポイントを推定する。聯立方程式を解いて求める。 この数列を既知のポイントとする。 こ…

TI-Nspire & Lua / 特殊な方程式に対する 7 段の Cowell 法

参考: パソコンで見る天体の動き, p.112-119 function cowell7(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] = ma…

TI-Nspire & Lua / 特殊な方程式に対する 4 段の Cowell 法 4 of 4 / 確かめる

参考: パソコンで見る天体の動き, p.112-119 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] = ma…

TI-Nspire & Lua / 特殊な方程式に対する 4 段の Cowell 法 3 / 前の近似値から次の近似値と近似値同士の最大誤差とを求める、それに意味のないとき

参考: パソコンで見る天体の動き, p.112-119 function approx2(funcs, w, preApprox, h) -- ({微分方程式}, {変化しない部分}, {前の近似値}, 刻み幅) local function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = math.abs(list…

TI-Nspire & Lua / 特殊な方程式に対する 4 段の Cowell 法 2 / 第 1 近似値、第 2 近似値、両者の最大誤差を求める、それに意味のないとき

参考: パソコンで見る天体の動き, p.112-119 function approx(funcs, w, temp, h) -- ({微分方程式}, {変化しない部分}, {假定値}, 刻み幅) local function maxOfErr(listA, listB) local sute = {} for i = 1, #listA do sute[i] = math.abs(listA[i] - lis…

TI-Nspire & Lua / 特殊な方程式に対する 4 段の Cowell 法 1 / 変化しない部分と假定値とを求める、それに意味のないとき

参考: パソコンで見る天体の動き, p.112-119 function staticVals(funcs, inits, h) --local function emptyMat(numOfRows) local sute = {} for r = 1, numOfRows do sute[r] = {} end return sute end local dim = #funcs local w = {} local f = {{},{},{…

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…