模造１（重力とエネルギー）
このページはマス旅の一部です。
今回は「重力とエネルギー」のシミュレーションをしよう。
geogebraのテキストでは上付きドット記号が使えないので（'）を使って時間微分の1階分を表すことにします。微分df/dtはf'で、d^2f/dt^2はf''とかきます。
１．解析解と数値計算解
課題１：「高さHtのビルの屋上から鉛直に初速度v0で玉を投げ上げる。
そのときの玉の運動を、横軸を時間ｔ、たて軸を位置yとして連続的に打点する。
重力加速度はg=9.81。」
運動方程式はF=ma=mv'=mx''（2階の微分方程式）だから、
自由落下の場合はF=mg=mv'=my''となるね。
＜解析解を利用する＞
玉の位置y、速度v、加速度は上向きを＋とします。
y,vはtの関数として解析的に書ける。
速度v=g(t)=v0-g*t
vグラフは上底v0,下底v0+g*t,高さtの台形となるので、面積は
(v0+v0-g*t)t*1/2=v0*t-1/2g*t^2
yの初期値Htとたすと、y=Ht+v0*t-1/2g*t^2
このyの式はHt,v0全部載せ
だから、これで曲線y=f(t)が完成。あとは、geogebraで作れるね。
＜オイラー法を利用する＞
時間ごとの変化に着目すると、
加速度が速度の瞬間増分（微分）だから、v'=g
速度が距離の瞬間増分（微分）だから、y'=v
時間増分th=0.1などとする。
初項v0, 公差 -g*th の等差数列になるから、リストvsができる。
v0, v0 -g*th, v0 -g*th*2,....
初項Ht, 階差vi*th の階差数列になるから、リストysができる。
Ht, Ht+v1*th, Ht+v1*th+v2*th, Ht+v1*th+v2*th+v3*th
.......
リストのi番目はHt+(v1+v2+...+v_{i-1})*th
＜geogebraのコード＞
g=9.81
th=0.1
rev_th=1/0.1
Ht=Slider(0, 10,1) #オブジェクトを画面に固定しましょう。
v0=Slider(-10, 10,1) #オブジェクトを画面に固定しましょう。
t=Slider(0, 5, 0.1) #設定＞上級＞アニメーション＞増加（51段階表示で連続表示）
#ここから解析解
f(x)=Ht + v0 * x-1/2 g * x^2
ex=If(f(t)>0,f(t),0)
P=(t,ex)) #地上で玉は止まる。
#ここからオイラー法
N=60*rev_th
vs=Sequence(v0 -g*th*k,k,0,N)
#階差部分を等差数列の部分和によってつくる。番号のズレ補正のため、初項を０にする。
vith = Join({0},Sequence(Sum(Take(vs,1,k)) *th , k,1,N-1))
ys=Sequence(Ht+Element(vith , k), k, 1,N)
idx = t * rev_th + 1 # (1..51）のレンジの離散量
ele=Element(ys,idx)
eu=If(ele>0,ele,0)
Q=(t, eu) #オイラー法の玉の動き
#Pの正確な動きと、Qの近似的な動きを比べてみよう。
thが0.1荒いので誤差が拡大してしまってますね。
thを細やかにすると点Pと点Qのズレは縮小するでしょう。
v0=0にすると自由落下になり、v0＜０にすると鉛直投げ下しになります。
鉛直投げ上げの解析解とオイラー法
２．NSolveODEを使おう
＜NSolveODEで課題１を解く＞
geogebraにはルンゲ＝クッタ法を内蔵した数値解を求める関数NSolveODEがあります。
それを使ってみよう。
距離の微分y'=v
速度の微分v'=-g(上向きが正）
解きたいの微分はy',v'で
そのy,vの初期値はy=Ht, v=v0です。
調べる範囲はさっきは指定しませんでしたが、12秒もいりません。t=10で十分です。
そこで、次の6行を入れます。
#定数
g=9.81
#初期値
Ht=Slider(0, 10,1)
v0=Slider(-10, 10,1)
#微分方程式
y'(t, y, v) = v
v'(t, y, v) = -g
#情報をまとめて投げるt=0からt=10まで調べる。t=0のときの初期値が`{ }の中。
NSolveODE({y', v'}, 0, {Ht, v0}, 10)
エンターすると、
一瞬で、y=f(t),v=g(t)の2関数を数値的に解明してくれます。
yのグラフを選んで、オブジェクト上の点を取ると、Aになります。
Aの設定を上級＞アニメーション＞速度0.3、表示条件y（A)、表示＞値にすると、
yのグラフにそって玉が地面まで移動して消えます。上級＞アニメーション＞増加にしないと、
点がもとに戻ってしまいますね。
鉛直投げ上げの微分方程式を数値的に解く
＜課題２：課題１に空気抵抗を追加する＞
空気抵抗があると加速度が-kvだけ変化して、a=-g-kvになるのです。
だから、微分方程式は、y'=vとv'=-g-kvの連立方程式になりますね。
空気抵抗のある落下運動をシミュレーションしようとすると大変なことになります。
それは、解析解ではv’変数分離法で積分でやってもvが指数関数になるので、
yは指数関数の積分になり、面倒な式になるからです。
オイラー法でもgeogebraで離散的にリスト処理でやろうとすると、
v自体が等差数列ではなくなるので、プログラミングが不可能ではないでしょうが、超面倒なのです。
もちろん、python,javascriptでwhile文を使えばできるでしょうが、大掛かりです。
でも、大丈夫！
ルンゲ＝クッタ法を内蔵した数値解を求める関数NSolveODEがgeogebraにはあります。
空気抵抗係数k=0.15という条件を追加して調べてみよう。
#定数
g=9.81
k=Slider(0,5,0.05) #追加最初はk=0.15にしましょうか。
#初期値
Ht=Slider(0, 10,1)
v0=Slider(-10, 10,1)
y'(t,y,v) = v
#この次の1行を変えるだけです。
v'(t,y,v) = -g - k*v
NSolveODE({y', v'}, 0, {Ht, v0}, 10)
k=0にすると、抵抗ゼロ、
kを上げていくと重力がキャンセルされて、
等速運動になるという
面白い現象がシミュレーションできます。
すごいぞ、geogebra!