はじめに

前の記事でメカニズムの基本を非常に簡単に検討したので、実践に移りましょう。頭で考えると簡単な理論で十分です。

したがって、タスク：

石は、h = 100 mの高さから初速度なしで垂直に投げられます。空気抵抗を無視すると、時間の経過に伴う地球の表面上の石の高さの関数として石の動きの法則が決まります。 10 m / s ^2に等しい重力加速度

簡単なタスクですか？はい、初歩的で、分析的ソリューションを備えており、読み書きのできる学生でも簡単に書くことができます。しかし、この単純なタスクは、私たちにとって非常に明らかな例となります。

1.タスクの形式化

シミュレーションの始まり。形式化により、進行中のプロセスを説明する数式の受信が理解されます。この場合、仮定が定式化されます。影響を無視できる要因によるモデルの簡略化のリスト。

仮定は、このタスクに適用可能であり、それに応じて、石は重要なポイントと見なすことができます。 1つの力がこのポイントに適用されます-重力なので、空気抵抗がないという仮定を使用します。

次の仮定は、私たちが石を投げる高さが惑星のサイズと比較して無視できるため、地球が平らであり、したがってその表面の曲率が無視できるためです。そして、重力は一定であると考えられ、地球の表面に垂直に向けられ、大きさが等しい

P = m \、 g

$P = m \、g$

ここで、gは地球の表面での重力加速度です。次に、石の運動方程式を作成します。これらの方程式を覚えておいてください

b e g i n c a s e s ＆ m \、 d d o t x = s u m F_{k x} ＆ m \、 d d o t y = s u m F_{k y} ＆ m \、 d d o t z = s u m F_{k z} \終 了 ケ ー ス

$\ begin {cases}＆m \、\ ddot x = \ sum F_ {kx} \\＆m \、\ ddot y = \ sum F_ {ky} \\＆m \、\ ddot z = \ sum F_ {kz} \終了{ケース}$

それらの左側の部分にはまだ興味がありませんが、右側には座標軸上の点に適用される力の投影の合計が含まれています。 x軸とy軸を表面上に配置し、z軸を表面に垂直に上向きにします。力は1つだけで、x軸とy軸上の投影はゼロに等しく、z軸上の投影は負です。これは、力が軸の方向に向かっているためです。

b e g i n c a s e s ＆ m \、 d d o t x = 0 ＆ m \、 d d o t y = 0 ＆ m \、 d d o t z = - m \、 g e n d c a s e s

$\ begin {cases}＆m \、\ ddot x = 0 \\＆m \、\ ddot y = 0 \\＆m \、\ ddot z = -m \、g \ end {cases}$

石の質量は明らかにゼロに等しくないため、結果の方程式の両側を安全に分割できます

b e g i n c a s e s ＆ d d o t x = 0 ＆ d d o t y = 0 ＆ d d o t z = - g e n d c a s e s

$\ begin {cases}＆\ ddot x = 0 \\＆\ ddot y = 0 \\＆\ ddot z = -g \ end {cases}$

私は気にしません、石の動きが厳密に垂直に発生することを証明していますが、これは正式な観点から行わなければなりません。最初の2つの方程式の右側のゼロは、これらの軸に沿った移動が不可能であることを意味しません。最初のニュートンの法則を思い出してください。これについては、次の記事で詳しく説明しますが、今のところ、最終微分方程式を書き出すことで、運動の1次元性を正しく設定します。

d d o t z = - g

$\ ddot z = -g$

私たちが成功したという事実は少なからずありません-問題で起こっていることのプロセスの数学的モデル。パトス？

いやこの方程式を分析すると、たとえば、この方程式には質量がないため、石の質量はその運動の法則に影響を与えないと結論付けます。方程式を解かなくても、私たちはすでにペンと真空中の鉛片を使った実験の有効性を公式に証明しており、彼らはそれを学校で見せたがっています（月にそれを繰り返した人もいます）。

分析ソリューションを取得するのは簡単です、私も気にしません、それはそうです

z （ t ） = h - f r a c g \、 t^{2} 2

$z（t）= h-\ frac {g \、t ^ 2} {2}$

しかし、これを数値的に解決する方法は？そして、それはどういう意味ですか-「数値的に」？

2.一階微分方程式の数値積分

最初の注文は何ですか？前回、運動方程式は2次であると言いました。すべて正しいのですが、コンピューターで拡散を解くためのほとんどの方法は、1次方程式しか解けません。 2次方程式を直接積分する方法（Verlet法など）がありますが、現在はそうではありません。

第一に、この方程式は、順序を減らすことができるようなタイプです。右側は未知の関数に依存しないため（zはありません）、そのことを思い出してください。

d d o t z = d o t v_{z}

$\ ddot z = \ dot v_z$

z軸の加速度の投影は、同じz軸の速度の投影の1次導関数に等しくなります。いいね

d o t v_{z} = - g

$\ dot v_z = -g$

これが一次方程式です。この番号は必ずしも通り抜けるわけではありません（コーシーのフォームについては今はお話ししません！）が、この場合はすべて正常です。座標ではなく、ポイントの速度を検索します。次は？そして

f r a c d v_{z} d t = - g

$\ frac {dv_z} {dt} = -g$

なぜなら、導関数は、関数の無限小増分（速度）と、それを引き起こした引数の無限小増分（時間）の比だからです。ほんの少しの時間しかかかりません

\デ ル タ t

$\デルタt$ 考慮できるほど小さい

f r a c d v_{z} d t \約 f r a c D e l t a v_{z} D e l t a t

$\ frac {dv_z} {dt} \約\ frac {\ Delta v_z} {\ Delta t}$

なに？そして、これは何ですか

D e l t a v_{z} = - g \、 D e l t a t

$\ Delta v_z = -g \、\ Delta t$

速度が上がりました。負の増分。このように、落下する石は同じように加速します！はい、そうです。その速度、その速度のベクトルは下向きになります。したがって、z軸上のこのベクトルの投影は負になります。そうです、絶対値で下向きのベクトルの投影を取得します。速度の初期値はゼロであり、つまり

v_{z}^{（ 0 ）} = 0

$v_z ^ {（0）} = 0$

速度の増分を計算できるという事実を利用して、0.1秒で速度がどのようになるかを計算します

v_{z}^{（ 1 ）} = v_{z}^{（ 0 ）} - g \、 D e l t a t = 0 - 10 c d o t 0.1 = - 1.0

$v_z ^ {（1）} = v_z ^ {（0）}-g \、\ Delta t = 0-10 \ cdot 0.1 = -1.0$

そしてさらに0.1秒後

v_{z}^{（ 2 ）} = v_{z}^{（ 1 ）} - g \、 D e l t a t = - 1.0 - 10 c d o t 0.1 = - 2.0

$v_z ^ {（2）} = v_z ^ {（1）}-g \、\ Delta t = -1.0-10 \ cdot 0.1 = -2.0$

そしてさらに0.1秒後

v_{z}^{（ 3 ）} = v_{z}^{（ 2 ）} - g \、 D e l t a t = - 2.0 - 10 c d o t 0.1 = - 3.0

$v_z ^ {（3）} = v_z ^ {（2）}-g \、\ Delta t = -2.0-10 \ cdot 0.1 = -3.0$

うーん、それでかなり長い間続けることができますが、1秒の時間間隔に制限します。

時間s	速度m / s
0.0	0.0
0.1	-1.0
0.2	-2.0
0.3	-3.0
0.4	-4.0
0.5	-5.0
0.6	-6.0
0.7	-7.0
0.8	-8.0
0.9	-9.0
1.0	-10.0

つまり、式を使用して

v_{z}^{（ i + 1 ）} = v_{z}^{（ i ）} - g \、 D e l t a t

$v_z ^ {（i + 1）} = v_z ^ {（i）}-g \、\ Delta t$

ポイントの速度の時間依存性が得られました。そして、あなたがする必要があるのは、現在の瞬間の速度値を取得し、そこに速度が新しい瞬間に受け取る増分を追加することです。

\デ ル タ t

$\デルタt$ 秒。ここでは、時間の増分を統合ステップと呼びます。そして、現在の時間に目的の関数の導関数の値にステップを乗じて増分を計算します。ただ？はい、もちろんです。そして、私が書いた式には、微分方程式の数値解法のための明示的なオイラー法という名前があります。これは、計算値の新しい値が以前の値に依存する場合の、いわゆる繰り返し式です。

しかし、地球上のポイントの高さはどうですか？はい、すべて同じです。

f r a c d z d t = v_{z}

$\ frac {dz} {dt} = v_z$

速度の投影は対応する座標の導関数だからです。速度はすでにわかっているため、この方程式にオイラー法の式を適用します

z^{（ i + 1 ）} = z^{（ i ）} + v_{z}^{（ i ）} \、 D e l t a t

$z ^ {（i + 1）} = z ^ {（i）} + v_z ^ {（i）} \、\ Delta t$

この式に従って、テーブルに別の列を追加します

b e g i n a l i g n ＆ z^{（ 1 ）} = z^{（ 0 ）} + v_{z}^{（ 0 ）} \、 D e l t a t = 100 + 0.0 c d o t 0.1 = 100 ＆ z^{（ 2 ）} = z^{（ 1 ）} + v_{z}^{（ 1 ）} \、 D e l t a t = 100 - 1.0 c d o t 0.1 = 99.9 ＆ z^{（ 3 ）} = z^{（ 2 ）} + v_{z}^{（ 2 ）} \、 D e l t a t = 99.9 - 2.0 c d o t 0.1 = 99.7 e n d a l i g n

$\ begin {align}＆z ^ {（1）} = z ^ {（0）} + v_z ^ {（0）} \、\ Delta t = 100 + 0.0 \ cdot 0.1 = 100 \\＆z ^ {（2 ）} = z ^ {（1）} + v_z ^ {（1）} \、\ Delta t = 100-1.0 \ cdot 0.1 = 99.9 \\＆z ^ {（3）} = z ^ {（2）} + v_z ^ {（2）} \、\ Delta t = 99.9-2.0 \ cdot 0.1 = 99.7 \\ \ end {align}$

などなど

時間s	速度m / s	高さ、m
0.0	0.0	100.0
0.1	-1.0	100.0
0.2	-2.0	99.9
0.3	-3.0	99.7
0.4	-4.0	99.4
0.5	-5.0	99.0
0.6	-6.0	98.5
0.7	-7.0	97.9
0.8	-8.0	97.2
0.9	-9.0	96.4
1.0	-10.0	95.5

うーん、まず第一に、速度は時間とともに変化するため、身長はすでに不均等に変化していることがわかります。現在、デリバティブ自体は時間依存です。しかし、すでに最初のステップで、何かがおかしいことに気づきました-速度はすでにそこにありますが、高さはまだ100メートルです。どうして？

これは、各ステップで微分（速度）定数を設定したために発生しました。このメソッドは、ステップ間でソリューションに何が起こるかについての情報を提供しません。したがって、誤差が累積し、得られた解と正確な解を比較します

z （ t ） = 100 - 5 \、 t^{2}

$z（t）= 100-5 \、t ^ 2$

時間s	速度m / s	高さ、m	正確な解、m
0.0	0.0	100.0	100.00
0.1	-1.0	100.0	99.95
0.2	-2.0	99.9	99.80
0.3	-3.0	99.7	99.55
0.4	-4.0	99.4	99.20
0.5	-5.0	99.0	98.75
0.6	-6.0	98.5	98.20
0.7	-7.0	97.9	97.55
0.8	-8.0	97.2	96.80
0.9	-9.0	96.4	95.95
1.0	-10.0	95.5	95.00

はい、私たちの石は空中で凍っているようです。数値解法は分析解法よりも遅れており、さらに考えると、アカウントのエラーが大きくなります。各ステップでますます大まかな近似を取るため、エラーが累積します。どうする

まず、ステップを減らすことができます。 10回言ってみましょう

\デ ル タ t = 0.01

$\デルタt = 0.01$ 秒

時間s	速度m / s	高さ、m	正確な解、m
0.0	0.0	100.0	100.00
0.1	-1.0	99.96	99.95
0.2	-2.0	99.81	99.80
0.3	-3.0	99.57	99.55
0.4	-4.0	99.22	99.20
0.5	-5.0	98.78	98.75
0.6	-6.0	98.23	98.20
0.7	-7.0	97.59	97.55
0.8	-8.0	96.84	96.80
0.9	-9.0	96.00	95.95
1.0	-10.0	95.05	95.00

カウント終了時のエラーは0.05メートルを超えておらず、これは以前の値の10分の1です。ステップをさらに10倍減らすと、さらに正確なソリューションが得られると想定できます。実際、このようなテーブルを取得するには、0.1刻みで10ポイントのみの値を表示して、10回ではなく100回の反復が必要です。ステップ0.001では、1000回の反復が必要になり、結果は

時間s	速度m / s	高さ、m	正確な解、m
0.0	0.0	100.0	100.00
0.1	-1.0	99.9505	99.95
0.2	-2.0	99.8010	99.80
0.3	-3.0	99.5515	99.55
0.4	-4.0	99.2020	99.20
0.5	-5.0	98.7525	98.75
0.6	-6.0	98.2030	98.20
0.7	-7.0	97.5535	97.55
0.8	-8.0	96.8040	96.80
0.9	-9.0	95.9545	95.95
1.0	-10.0	95.0050	95.00

これらの計算を手動で実行しようとした場合、高精度が必要な場合にどれだけ単調で時間がかかるかを理解できます。それが、数値シミュレーションの全盛期がコンピューターの出現と一致した理由です。これらは、数値に対して多くの単調な操作をすばやく実行するために必要です。

オイラー法は、微分方程式を積分するための最も簡単な既知の方法です。この単純な例から、メソッドのエラーは積分のステップに正比例することがわかります。これは実際にそうです。このようなメソッドは、1次精度メソッドと呼ばれます。

たとえば、4次精度のルンゲクッタ法を適用すると、ステップ0.1でも計算の精度が向上します。しかし、これは別の話です。

おわりに

それについて考えてみてください...私たちは非常に単純な例を見てみました。私たちはコンピューターさえ使用しませんでしたが、世界の非常に強力なスーパーコンピューターが宇宙の生活の初期段階をシミュレートする仕組みをすでに理解しています。もちろん、すべてがより複雑ですが、原則は同じです。

あなたが自分の手に渡る強力なツールを想像してください。これは、コンピューターを使用しない最後の記事です。オクターブを約束しました。次回は彼になります。

動的システムモデリング：ODEを解くための数値的手法

はじめに

1.タスクの形式化

2.一階微分方程式の数値積分

おわりに

More articles: