11回目　微分方程式の数値演算(1)

2015年度　中間試験の解説

シミュレーションの考え方
計算機シミュレーションは大変便利で，現在の様々な製品を開発する上でなくてはならないものに
なったと言えます．機械工学をバックグラウンドとする技術者にも大変有用なものであり，例えば，
3D CADに様々なシミュレーション機能を組み合わせて解析し，設計に役立てることができます．

例

• 機構解析と干渉（接触）
• 構造解析や振動解析
• マスプロパティと動力学
• 伝熱解析
• 流体解析

もはや機械を設計する上で実務上なくてはならない存在ともいえるのがシミュレーションであり，みな
さんの想像の通り，多くのパラメータ，数多くの要素を一度に演算できるのが大きなメリットです．
設計の初期段階では，仕様を決め，それを満たすための代表値を中心とした計算を中心に取り組む
わけですが，あまりにも多くのパラメータを数学的に綺麗な形で算出することは事実上困難なのです．
ところで，以下のようなキーワードを既に聞いていると思います．

• 解析的解法
• 数値的解法

解析的とは，数学的に綺麗，即ち，数式をきちんと立て，その数式を変形すると，解についての方程式
を示せるところまで数式に従って解ける方法と言えます．例えば，パラメータx0〜x9まであってyを知り
たいとき，y=f(x0〜x9)のように表現でき，ｙを得ることができるような場合を指していると言ってよいと
思います．

これに対し，数値的とは何かを説明します．
数学的に解くことには変わりがないのですが，ここで例えば，材料力学で扱うFEMを考えてみましょう．
材料は，連続体として見なせます．
細かく見ていき，組織あるいは分子まで分解し，それらの相互の関係を計算するという方法もあります
が，常の機械設計においてはそこまではやりません．結果に対する計算コストがかかり過ぎて意味が
ないからです．
単純な構造であれば解析的な解法が使えますが，複雑な形状における応力分布，とりわけ応力集中
などを知りたい場合があります．構造材をメッシュとして捉え，多数の点に関する方程式を立てるわけ
ですが，数千点以上になると，とてもそれらを一意に解くことなど出来ないことは想像に難くないと思い
ます．
詳細は割愛しますが，ある初期状態，ある条件を与えたとき，モデルを構成する多数のパラメータに関し，
一つのパラメータに関して仮に値を与えて解き，その影響をさらに計算して次のパラメータを解き，それを
波及させて解いていくことを繰り返したりします．
コンピュータは知っての通り離散量を扱います．すなわち有限個の要素やパラメータを用いていることが
特徴なのです．以前お話ししたサンプリングを思い出してみて下さい．時間軸で区切って標本化し，振幅
で区切って量子化しています．もともと連続量だったものを時間で区切り，その時の値も有限の段階で
示します．つまりは，あらゆるパラメータに関して，段階的にしか扱うことができず，FEMのように少しずつ
変化させてその様子を観察するという方法を用いるのです．
このような性質から，多数の観測点について連立した方程式を多数解くために有効なものと考えて利用
します．

前置きが長くなりましたが，数値的であることの問題は，中間試験の最後の問題でも答えていただいた
通りなのです．以下にその説明を書いていきます．

切り捨て，丸め誤差
皆さんが使っている変数の型floatやdoubleは約15〜17桁位の精度を持っています．
それ以上の値を入れても勝手に切り捨てられたり，丸められたりします．
数学・物理的には同じ定数であっても，小数点表記だけで表記するのか，10の冪乗と分けて表記するか
によっても，計算精度に違いが生じます．

誤差の蓄積
計算機のお得意は単純な計算を早く繰り返して行うことです．
シミュレーションはその性質を利用して，値を少しずつ変化させて，方程式を解くことを繰り返します．
動力学を解くための連立微分方程式を解く方法など代表例かと思われます．
極端に小さい値と極端に大きな値を掛けたり割ったりすると誤差が生じます．
切り捨てや丸め誤差に起因してしまいます．
またある程度注意深く数値を選択しても，誤差を含んだ計算を繰り返すうち，それが蓄積してしまうことも
あるため，さらに大きな誤差となり，シミュレーションの結果の妥当性が危ぶまれる結果につながります．

シミュレーションの条件
シミュレーションのねらいは，複雑な計算に対し正確な解を早く得られることに他ならないと思います．
計算機負荷に関しては従来は負担になるという都合もあり，刻み幅を大きくしたり，メッシュを大きくしたり
するなどして，データ量や計算量を減らしてきたが，今となってはむしろ，不具合の起きない条件を探す
必要があります．計算のステップを粗くしたところで誤差を発生することもありますが，細かくし過ぎても
やはり解がおかしくなることがあるのです．

参考
シミュレーションの例を挙げると，微分方程式を解くための数値演算を行う場合があります．
時間刻みごとの解の組を計算するわけですが，時間刻みを小さくすることで計算量が増えることを避ける
ため，刻みを粗くすると，困ったことに解が発散するという問題が存在してしまいます．
もちろん刻みを小さくすればきちんと解を持つので，微分方程式の解がないわけではないため，問題は，
刻みを粗くしすぎたことと考えられます．
時間刻みを闇雲に小さくすれば，精度の細かい計算だから良いかというと必ずしもそうではありません．
時間刻みを極端に小さくすると変化が僅少になるため，切り捨てられてしまった故の誤差が発生してしまい
これを繰り返すと誤差を累積することとなるのまた別の問題として発生してしまうのです．
ですので，シミュレーション条件を変えたりする必要があるとともに，解析的に解いた場合の大まかな傾向を
知っておかないとシミュレーションが正しいかどうかも分からなくなってしまうことがあります．

参考）　PCやCPU（個人ユーザー向けのちょっとハイスペック）の変遷を紹介します．

大雑把に調べたものですが，様々なアーキテクチャ（CPUの構造）が開発されクロックがどんどん
上がってきたた上に，近年はコア数が増え，並列処理的な処理がなされるようになってきました．
内部の処理，周囲に配置されたキャッシュメモリなどの改良によって，単純なクロック数以上の
進歩を続けています．
10数年以上昔は，数値シミュレーション，CGなどは一晩，あるいは数日かけて計算させていたもの
ですが，いまやそんなこともなくなっています．
古い文献などは，ヘリコプタの制御ができる処理能力はとても実現できないと書かれていましたが，
今となってはRT（実時間）をはるかに超える処理能力を持ったものとなっています．

以上のCPUの進歩とメモリの大容量化（30年前は1Mbyte以下）を考えるとプログラムの処理能力
については考え方が変わった面もありますが，組み込み機器（制御機器や家電，設備を動かす
CPU）については，省電力とコストの問題から，CPUも数十MHｚオーダのものであり，メモリも大きく
ても十数Mbyte程度ですので，プログラムの効率化は今でも大変重要なのです．

常微分方程式の数値解法
今回は，一般的に使われながら，より精度の高いルンゲクッタ法と呼ばれるもので，1階の
常微分方程式dx(t)/dt = f(t, x)を解く際に利用できます．
時間刻みΔtを便宜上hとおくとして，h後の時刻と解は次のようになります．

ti+1 = ti + h
xi+1 = xi + k
k=(1/6)(k0+2k1+2k2+k3)

最後の式がルンゲクッタの解法の特徴ですが，kはそれぞれ，

k0 = h f(ti, xi)
k1 = h f(ti + h/2, xi + k0/2)
k2 = h f(ti + h/2, xi + k1/2)
k3 = h f(ti + h, xi + k2)

となります．複雑な計算をしているように見えますが，離散的ながら，重みを付けて4次の曲線として
数値解を出していると考えていただければよいかと思います．
なお，アルゴリズムは，

1. 初期値t0, x0を決定し，時間刻みhを適宜決めます．
2. k0 = h f(ti, xi)を計算します．
3. k1 = h f(ti + h/2, xi + k0/2)を計算します．
4. k2 = h f(ti + h/2, xi + k1/2)を計算します．
5. k3 = h f(ti + h, xi + k2)を計算します．
6. k=(1/6)(k0+2k1+2k2+k3)を計算し，次のステップ ti+1 = ti + hのkとします．
7. xi+1 = xi + kでh刻み後の解が算出できます．

あとは2〜7を繰り返します．for文などで繰り返せばよいかと思います．
以下にソースファイルがありますので，適宜利用して下さい．

ルンゲクッタ法→rkm.c
ルンゲクッタ・ジル法→rkgm.c
ミルンの方法→mlnm.c
アダムス法→admsm.c

いずれもルンゲクッタのアルゴリズムを基本としていますが，係数を変えたり，計算を追加して，精度を
上げる工夫をしたものです．但し，計算量は少しとはいえ増える傾向にあります．

2階微分方程式

ルンゲクッタ法は1階の微分方程式を解くための方法です．工夫すれば，連立した1階の微分方程式を
解くことができます．
機械工学で利用する微分方程式は2階もしくはそれ以上の高次のものが大半です．これを解くにはどの
ようにしたらよいかというと，この連立して解く方法を利用するのです．
dx/dt=zとおいて，導関数を新たに変数にしてしまうことにします．x'= zだからx"= z'となるわけです．
従って，新たに変数が発生したため，

dx/dt = f(t, x, z)
dz/dt = g(t, x, z)

の連立方程式を立てて同時に解けばいいのです．

ti+1 = ti + h
xi+1 = xi + k
zi+1 = xi + l
k=(1/6)(k0+2k1+2k2+k3)
l=(1/6)(l0+2l1+2l2+l3)

但し

k0 = h f(ti, xi, zi)
l0 = h g(ti, xi, zi)
k1 = h f(ti + h/2, xi + k0/2, zi + l0/2)
l1 = h g(ti + h/2, xi + k0/2, zi + l0/2)
k2 = h f(ti + h/2, xi + k1/2, zi + l1/2)
l2 = h g(ti + h/2, xi + k1/2, zi + l1/2)
k3 = h f(ti + h, xi + k2, zi + l2)
l3 = h g(ti + h, xi + k2, zi + l2)

のような計算をすればいいのです．式は増えていますがやっていることは変わりません．

補足説明
ミルンの方法とアダムスの方法について，補足説明をします．
ルンゲクッタ法は基本になりますが，4次までの近似で，なおかつ多角形的なものであることはオイラー法と
変わりがありません．
その欠点を経験則的に解決したのがルンゲクッタ・ジル法であり，係数が細かく調整してあります．それでも，
微分方程式の曲線のカーブの変化が大きい時などは誤差が蓄積します．
ミルンの方法，アダムスの方法はこの問題点を解決する方法であり，関数中に設けたepsというパラメータが
あることに着目してみて下さい．
これは設定した誤差であり，一ステップ当たりの誤差がeps以下になるまでルンゲクッタによる解を修正すると
いうものです．
ルンゲクッタによる４ステップ分の計算をまず行い，これらの解を基に修正を１ステップごとに繰り返していくと
いうものです．ですから計算量は多くなります．epsを小さくすればするほど計算量が増える可能性が大きく
なるといえます．
ミルンの方法→comp_mlnm.c
アダムスの方法→comp_admsm.c

数値演算については精度や型のことを知っておく必要があります．以下のページを必ず読んでおいて下さい．

変数の精度について

課題１　以下の関数を作成しし，上記の4つの方法によって常微分方程式を解け．
関数に関するヒントはここにある．
4つのファイルをダウンロードし，ソースファイルを完成せよ．また解いた結果を比較し，考察せよ．

x(t)' = - b * x(t) * x(t) - a * x(t) + g

これは速度に比例する抵抗と速度の2乗の抵抗を受ける物体が落下する系に関する運動方程式である．
定数はa=0.1 [1/sec], b=0.02 [1/m], g=9.807 [m/sec/sec]とするがa,bは任意．
t=0〜任意の時間について解けばよい．

xの時間関数であるから以下のような引数として，関数を作ればよい．

double fnc(double t, double x)

課題２　以下の関数について，ルンゲクッタ法を用いて常微分方程式を解け．

x"(t) = - (c/m) x'(t) - (k/m) x(t) + u/m

ばね-マス-ダンパー系でc/m = 0.02[N sec/m/kg], k/m = 0.1[N/m/kg], u/m = 1[N/kg]位をパラメータとし，
初期値はx'(t)= 1[m/sec], x(t) = 0.05[m]として計算してみよ．
数値は変えてもかまわない．

ルンゲクッタ法→rkm.c　をダウンロードし，名前を変えて先のアルゴリズムに合わせて変更すればよい．
要は，上記のアルゴリズムをそのまま並べて書けば充分である．rkm_2.cのような名前にしておけばよい．
運動方程式の関数は dx/dt = f(t, x, z)，dz/dt = g(t, x, z)となることに注意せよ．

課題に取り組むにあたり，ソースファイルは以下のものを参考にしてもよい．
ミルンの方法　→comp_mlnm.c
アダムスの方法　→comp_admsm.c

補足説明

提出物
・ソースファイル
・グラフを作成したファイル（課題1と2ともt-xの関係が分かるように記すこと。）
・考察を書いたテキストまたはワードファイル

シミュレーション時間は以下のようにするとわかりやすいグラフとなるのでこの範囲で作成すること。
課題１　20秒程度
課題２　50〜100秒程度

2階の微分方程式を解くヒント

元の運動方程式について
m x" + c x' + kx = u
z = x'
とおくと、
z' = x"となるため、z'について解くと
以下の連立方程式となる。

dx/dt = x' = z
dz/dt = z' = -(c/m)z - (k/m) x + (u/m)

これらを関数として代入すればよい。