微分方程式を解く

ローレンツアトラクタの節の最初に例として示したMathematicaとGrapherの図では、微分方程式を入力してそれを解かせています。MathematicaではNDSolveという微分方程式を解く関数を用いています。Grapherでは微分方程式と条件（初期値）を書くだけでグラフを表示することができます。Grapherにはオプションがあって、「オイラー法」と「ルンゲクッタ法」のいずれかが選べるようになっています。Mathematicaでは、NDSolveのチュートリアルに、次のように書かれています。

「NDSolveでは解を求めるために反復手法が使われる。それは、xのある値を始点とし、一連の刻みをとりながらxminからxmaxの範囲をくまなく探索していく」

この、「反復手法」に、オイラー法とルンゲクッタ法があるのです。ここでは、オイラー法を用い、ルンゲクッタ法に詳しく言及することはしません。ルンゲクッタ法については、ウィキペディアなどを参照してください。高校数学の教科書にある「微分」「微分方程式の解法」から始めて、Cinderella.2で実験しながら、カオスを探究するのに必要な方法を考えていきましょう。

微分の定義と記法

まず、微分の定義と記法、その意味をグラフで復習しておきます。

関数(x)の点(a,f(a)) における微分係数を次のように定義します。

このaを変数xにすれば、そのまま導関数の定義になります。

ここで、前者の記法 f'(x) エフダッシュ（イギリスではダッシュではなく「プライム」という）はラグランジュが用いた記法、後者の dy/dx はライプニッツが用いた記法です。ここで、ライプニッツの記号が分数形になっていることに注意しましょう。また、ニュートンはｙの上に点（ドット）を打って表しました。

微分係数は、図形的には、曲線 y=f(x) の点(a,f(a))における接線の係数となります。したがって接線の方程式は y=f'(a)(x-a)+f(a)です。

そこで、(a,f(a))における接線を短い線分として表示するスクリプトを作ってみましょう。点(a,f(a))を自由に動かせるようにするため、文字も変数らしい t にします。

まず準備として、座標軸と方眼を表示しておき、「直線を加える」モードでx軸上に直線ABを引き、「点を加える」モードで直線AB上に点Cを取ります。「要素を動かす」モードで、点Ｃが直線ABすなわちx軸上を動くことを確かめておきましょう。

準備ができたら次のスクリプトを書きます。

ーーーーーーーーーーーーーーーーーーーーーーーーーー

f(x):=x^2;

f1(x):=d(f(#),x);

s(x,t):=f1(t)*(x-t)+f(t);

t=C.x;

plot(f(x));

plot(s(x,t),start->t,stop->t+0.5,color->[1,0,0]);

ーーーーーーーーーーーーーーーーーーーーーーーーーー

１行目は関数の定義。ここでは簡単な2次関数です。

２行目は導関数の定義。ただし、limをとる関数はないので、CindyScriptにある微分の関数 d() をそのまま使います。

３行目は接線の方程式の定義

１行空けて、4行目では点Cのx座標を変数tとして取り込んでいます。

５行目は関数f(x)のグラフを描画

６行目で点(t,f(t))から長さ0.5の接線を赤で表示します。

点Ｃを動かすと、接線も動きます。これで、「点の近くでは、曲線と接線はほぼ一致する」ことを体感しましょう。

ということは、接線の傾き＝微分係数を使って、曲線上である点に近い点すなわち関数の値が求められることがわかります。

こうして、すぐ近くの点を次々にとっていけば、導関数からもとの関数の値が求められることになります。

これが、原理です。

では、この原理に基づいて関数値を求めていったものと、関数の実際の値を比較するスクリプトを書いてみましょう。

導関数から関数値を求める

スクリプト例です。

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

f(x):=1/2*x^2;

f1(x):=x;

dx=0.1;

g(x,y):=y+dx*f1(x);

x=0;

y=0;

plot(f(x));

repeat(40,

draw([x,y]);

y=g(x,y);

x=x+dx;

);

ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー

１行目、もとの関数の定義です。係数が1/2 の２次関数。

２行目、その導関数です。

３行目、ｘの増加量です。

４行目、導関数を用いて、あるx,y から次のyの値を計算します。この式の意味が問題です。導関数の値すなわち微分係数は接線の傾きでした。これを用いて、点(x,y)に近い (x+dx , y ) の yの値を近似的に求めています。あとで少し異なる切り口で確かめますが、なぜこの式で近くの関数値が近似的に求められるのか考えてください。

５，６行目出発点の x , y の値です。初期値といいます。

７行目、関数のグラフを描きます。これが「本物」の値のグラフです。

８行目〜１２行目、すぐ近くの値を計算する、という操作を40回繰り返します。ｘの増加量が0.1ですので、x=4までの値が0.1刻みで求められることになります。

９行目、計算した点をプロットします。

10行目、xが0.1だけ増えた点のyの値を計算します。

11行目、xの値を0.1だけ増やします。

実行結果は次のようになります。緑の点が、導関数から計算したもとの関数の値を示す点です。

スクリプトの１行目の関数の定義、２行目の導関数（自分で計算して定義するのですよ）、３行目のxの増加量、５，６行目の初期値をいろいろ変えて実験してみましょう。

微分方程式から差分方程式を作る

導関数を表すのに、ここではライプニッツの記号が便利ですのでこれを使います。

先ほどの例は、

から

の関数値を計算しました。

①の式は、導関数の式ですが、これを方程式と考えて微分方程式といいます。ここから、もとの関数を求めることを「微分方程式を解く」といいます。この例はもっとも簡単なもので、積分すれば①から②は求められます。そのときに、積分定数Ｃが出てきますが、「x=0のときy=0」という条件をつければ、この積分定数の値は決まります。この条件を初期条件といいます。

さて、先ほどやったのは何かというと、積分をしないでもとの関数の値を求めた、ということです。積分を行なって求めるのを「解析的に解く」といい、値を次々に求めていってもとの関数がどんな関数かを調べることを「数値的に解く」といいます。先ほどの例でわかるように、数値的に解いた場合、解析的に解いた「解」と、完全には一致しません。近似的に解けた、ということです。

では、①の右辺がyの場合の微分方程式を解いてみましょう。高校数学では発展的内容となっています。初期条件を x=0のときy=1とすると