物理の小さい余白 物理の大事な骨の部分だけを独自にまとめましたサイトです
ご意見ご感想はこちらまで（「」内に@gmail.comを付けてください）　「ysdswallow」
top page 微分・積分のまとめ ニュートン力学のまとめ
電磁気のまとめ 解析力学のまとめ 数値計算（微分方程式）のまとめ
おすすめ参考書等紹介

[次へ][目次]

Chapter 1
差分法 (Euler 法)

1.1 Euler 法

物理の基本的な式はみな微分方程式です。ニュートンの運動方程式も量子力学のシュレディンガー方程式もそうです。そ れで、ここではコンピュータを使ってどうやって微分方程式を計算するのかを簡単に紹介したいと思いま す。

さて、最も簡単な形としてオイラー（Euler）法というものがあります。まずもっとも簡単な微分方程式として例え ば

dy-= y dx

のように y の x での微分が x,y の関数でかけているとき

dy- dx = f(x,y)

を考えていきます。このままでは当然コンピュータには計算させられません。そのためにはなんとか電卓で計算できる形、 つまり四則演算の形にしてやらなければいけません。さて、それにはなにより微分のところがどうにもなりませんね。そも そも微分とは

dy-= y(x+-Δx-)--y(x-),Δx → 0 dx Δx

でした。これから０にはできないけども出来うる限り小さいΔx を使えば

y(x-+-Δx-)--y(x)- Δx = f(x,y) (1.1) y(x + Δx) = y(x) + f(x,y)Δx (1.2)
となりこの計算を n 回繰り返していけばy(x+ nΔx ) といったように y をどんどん求めていくことができます。これが 「Euler 法」です。実に簡単ですね。ためしに
dy- dx = y,y(0) = 1

で計算してみましょう。この式の解は

x y = e

なので、今x = 0 から１まで計算するとなると本当の答えはy(1) = e = 2.71828 となるはずです。さて、この計算のた めの C 言語でのプログラムのソースはこちらです。

/*Euler法の練習*/  
 
#include <stdio.h>  
 
int main(void){  
 double y=1.0; /*yの初期値は１です*/  
 double x=1000; /*この逆数がメッシュの大きさです*/  
 long i;  
 
 for(i=1;i <= x;i++){  
  y += (1/x) * y;  
 }  
 printf("%f\n",y); /*結果を表示します*/  
 return 0;  
}

このソースでは基本的に

yn+1 = yn + ynΔx,(yn = y(nΔx ))

を繰り返しです。Δx = 0.001 にして、 x=0 から１までなので上の式を 1,000 回繰り返したのですが、ここで問題なの がΔx の値（以下メッシュという）をいくつにするかです。理論的にはメッシュをできるだけ小さくしたほうが 実際の微分に近づくわけだからそのほうがいいんですが、当然０にはできないしまた、メッシュを小さく すればするほど計算回数が多くなり時間がかかってしまいます。Δx = 0.001 のこのソースでは結果は 2.716924 と誤差は 0.05 ％というところです。メッシュの大きさと真値との誤差をグラフにしたのが次の図で す。

PIC

Figure 1.1: メッシュの大きさと誤差

です。縦軸に誤差の％の対数、横軸がメッシュの幅の対数です。これからともに対数を取ったものは直線になっ ていることがわかると思います。この場合はかなり幅の小さいメッシュにしても時間はそれほどかかりま せんが、もっと複雑な式の場合は誤差と時間の兼ね合いを考えてメッシュの幅を決める必要があるでしょ う。

1.2 二次の微分方程式

さて、代表的な物理の方程式といえばニュートンの運動方程式

d2r F = m dt2

ですが、これは二回微分の式です。これに限らず二回微分の方程式は物理では大変多いです。一般に二回微分の 式

ｨy = f(y,x)

を Euler 法で計算するにはまず、ｙの一回微分を新たにv(y,x) = ˙y という関数だと考え

˙v = f(y,x) ˙y = v(y,x)
をいうふたつの一回微分の式に分解してそれぞれ計算していけば、求めることができます。この方法を使えば基本的に何次 元の微分方程式でも計算可能です。さて試しにオイラー法で運動方程式を解いてみましょう。計算する式とし て重力による惑星の軌道を考えてみましょう。簡単のため重力定数や質量など仮に１だとした式を考える と

d2x- x- dt2 = - r3 d2y y dt2-= - r3
となります。ただし
x(0) = 1,y(0) = 0 dx(0) dy(0) -dt--= 0,-dt--= 1 ∘ -2---2- r = x + y
です。さて、この解は楕円軌道だと仮定し
x = a cos ωt+ C1,y = b sinωt +C2

を代入して計算してやればすぐに半径１、中心０の円であることがわかります。ではさっそくこの式をオイラー法で計算し ましょう。ソースは変数が増えるのでややこしくなりますがこうです。

/*Euler法 メッシュの幅は０．２としｔ＝２０まで計算する*/  
 
#include <stdio.h>  
#include <math.h>  
 
int main(void){  
 double x =1;  
 double u[2] ={0};  
 double y =0;  
 double v[2] ={1}; /*u,vはそれぞれx,yの一回微分*/  
 int k;  
 
 for(k=0;k < 1000; k++){  
  u[1] =u[0]+ -0.02 * x / sqrt(y*y+x*x);  
  v[1] =v[0]+ -0.02 * y / sqrt(y*y+x*x);  
/*u,vは二回使うので一時的な保存先としてu(v)[1]を使う*/  
  x += 0.02 * u[0];  
  y += 0.02 * v[0];  
  u[0]=u[1];  
  v[0]=v[1];  
 
  if(k%20==0) printf("%f %f\n",x,y); /*20回に一回だけ結果を表示させる*/  
 }  
 return 0;  
}

このソースでの計算は

xn- vn+1 = vn - r3 Δt yn- un+1 = un - r3Δt xn+1 = xn + vnΔt yn+1 = yn + unΔt
を繰り返すだけです。メッシュの幅は０．２、ｔ＝２０まで計算しました。さて、この結果は次のようになりま す。

PIC

Figure 1.2: Euler 法による惑星運動

赤い線が理論値で緑が計算結果です。これははっきりいって悲惨ですね。もちろんメッシュを小さくすればマシにな りますがだんだん外に膨らんでいくのは変わりません。これはオイラー法の欠陥なのです。これで地球の 100 年後を計算すれば太陽系から離脱してしまいます。そこで次にもっと精度の高い計算方法を紹介しま す。

[次へ][目次]