常微分方程式の数値解法

Nspired Pythonによる単振り子のシミュレーション / アニメ化する

tns_files/pend_single_tns at master · ti-nspire/tns_files · GitHub import RK4_Class import time import ti_system as tsys import ti_draw as td from math import sin as SIN from math import cos as COS from math import pi as PI class PEND_SING…

Nspired Pythonによる単振り子のシミュレーション / 常微分方程式の解法

参考: パソコンで見る天体の動き, pp.26-29 まず微分方程式が解けるようにする。何もインポートしていないので普通のPythonでもそのまま走る。 # Python による常微分方程式の数値解法_改良版 / 古典的 Runge-Kutta 法 class RK4: def __init__(self, funcs,…

常微分方程式をVIリファレンスとして入力する方法

ブロックダイアグラムに[スタティックVIリファレンス]函数を配置する。 配置したそのアイコンを右クリックし、[パスを参照...]を選択し、目的のVIを選択する。 再度そのアイコンを右クリックし、[タイプ指定されたVIリファレンス]に✔マークをつける。 これで…

[ODE4次ルンゲ・クッタ法VI] (ODE Runge Kutta 4th Order.vi)

RK4_Ex.vi - Google ドライブ [X(変数名)]が従属変数名排列。 [X0]が従属変数初期値排列。 [時間]が独立変数名。デフォルトは"t"。 ここでは従属変数に{"x", "y"}を指定したので、解は左からx, yの順に並ぶ。 詳細ヘルプに「ルンゲ・クッタ法を使用して、初…

Python / 常微分方程式の数値解法

Python で常微分方程式の数値解を求める。聯立微分方程式であっても方程式ごとに個別に def で定義する。リアルタイムシミュレーションが目的であるので 1 刻みだけ計算する (1 刻み計算して描画、1 刻み計算して描画、...... を繰り返す)。 関連: TI-Nspire…

Python / 常微分方程式の数値解法 / 補外法 / 應用例 / 単振り子 / 時間領域のグラフを描く

関連: ログイン - はてな from myODE.Extrapo import * # myODE.Extrapo は補外法で 1 階常微分方程式を解くためのプログラム from math import sin as Sin from math import pi as Pi from matplotlib import pyplot as plt class Pend: def __init__(self,…

Python / 補外法 / 補外にかける点列を中点法で求める

関連: TI-Nspire & Lua / 補外法 5 / 補外にかける点列を中点法で求める - def midPoints(funcs, t0, inits, H, S): dim = len(funcs) RANGE = 2**S+2 listRange = [None for i in range(RANGE)] temp = [None for i in range(dim)] listS = [None for i in …

Python / 補外法 / Neville の算法で補外値を求める

関連: TI-Nspire & Lua / 補外法 4 / 何かを偶数次多項式で補外する / Neville の算法で補外値だけを求める / 刻み幅を順次 1/2 にした場合 / 複数の点列に適用する - def neville(listOfLists): dim = len(listOfLists) S = len(listOfLists[0]) output = […

Python / 古典的 Runge-Kutta 法 / 1 ステップだけ計算してみる

def rk(func, t0, init, h): f1 = func(t0 , init ) f2 = func(t0 + h / 2, init + h * f1 / 2) f3 = func(t0 + h / 2, init + h * f2 / 2) f4 = func(t0 + h , init + h * f3 ) return init + (h / 6) * (f1 + 2 * (f2 + f3) + f4) # 次の微分方程式を 1/8 …

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

参考: パソコンで見る天体の動き, pp.160-166 Excelによる数値解析―オイラー法のうまい使い方, pp.94-105, pp.145-158 L4 を中心に、無視できるほど軽い質点をいくつか配置する。その動きを補外法で計算する。 Extrapo = class() function Extrapo:init(func…

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 法

参考: パソコンで見る天体の動き, pp.162-166 Excelで操る! ここまでできる科学技術計算, pp.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 / 補外法 / 質点同士の接近時のみ正則化する / それに意味のないとき

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

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

参考: パソコンで見る天体の動き, pp.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で学ぶ微分方程式, pp.204-209テキストは古典的 Runge-Kutta 法を使っているが、ここでは、特殊な方程式に対する Nyström 法を使ってみる。 縦の点線が平衡点: CoupledVib = class() function CoupledVib:init(funcs, t0, inits, initsDot, h, nu…

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 法

参考: 工業力学 (機械工学基礎講座), pp.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.h…

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で学ぶ微分方程式, pp.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.fun…

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

関連: TI-Nspire & Lua / 常微分方程式の数値解法のまとめ (『パソコンで見る天体の動き』) / クロージャを利用する - 同じ処理を今度はクラス化する。 1. 古典的 Runge-Kutta 法 構文: rkClassic({funcs}, t0, {inits}, h [, numOfDiv])▶ソースコードを表示…

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

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

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

参考: パソコンで見る天体の動き, pp.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 / まとめ

参考: パソコンで見る天体の動き, pp.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 one…

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

参考: パソコンで見る天体の動き, pp.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, …