第31回 数値積分
Cinderella.2では、Mathematicaのような記号処理としての積分はできませんが、定積分の値を求めることはできます。数値積分とは具体的な値について定積分を計算することで、リーマン積分の原理に基づいた方法のほか、台形公式、シンプソンの公式を用いる方法などがあります。
リーマン積分と区分求積法
この方法は定積分の定義といってもよいでしょう。区間[a,b]をn個の区間に分け（等分でなくてもよい）矩形の面積の総和を考えます。矩形の高さを、分割した各区間の最大値で取り、n→∞としたときの極限を上積分、最小値でとったときを下積分といい、それが一致したとき、定積分の値とします。厳密な定義としてはもう少し条件が必要ですが、イメージとしてはこのようなものです。
高校の数学Ⅲの教科書に載っている「区分求積法」は、n等分にして、各区間の端点に対する値を使って矩形の面積の総和の極限を考えます。ここで、「面積」と書いていますが、実際には、区間の幅×関数値ですので、関数値が負ならば通常の意味での面積にはなりません。むしろ、正負はともかく数値計算の結果として定積分を定義し、その区間で常に関数値が0以上ならば面積と同等である、と考えるべきでしょう。
これらは積分の意味としては分かりやすいのですが、実際に計算するには収束が遅いのです。
収束が遅いというのは、分点をどんどん増やしていったとき、積分値に近づく度合いがゆっくりだということです。
そこで、数値計算で積分値を求めるためにいくつか工夫します。
台形公式
矩形の代わりに台形を使えば曲線の下の図形により近くなります。
上の図と見比べれば一目瞭然ですね。
では、これをCinderella.2でプログラミングしてみましょう。
まず、いくつか点を作図します。シミュレーション形式にするために、分点の数をスライダで変えられるようにし、また、定義域も変更できるようにします。次が、作図した要素の一覧です。
非表示になっているAA,BB,C,Dは、座標軸を描くためのものです。S,Tはスライダの両端、Pはスライダ上の点、A,Bは定義域の両端です。
こうしておいて、次のスクリプトをDrawスロットに書いて実行すれば、シミュレーションができます。
-----------------------------------------------------------------------------------------------
f(x):=x^3-2*x^2+2;
A.y=0;
B.y=0; // 点A,Bがｘ軸から離れないようにする。
sx=A.x; //定義域左端
ex=B.x; //定義域右端
n=round((P.x-S.x)*2); // スライダから分点の数を得る
plot(f(#));
dx=(ex-sx)/n; //刻み幅
xn=apply(0..n,sx+#*dx); // sxからexまでをn等分した分点のリスト
yn=apply(xn,f(#));
sum=0;
repeat(n,
drawpoly([[xn_#,0],[xn_#,yn_#],[xn_(#+1),yn_(#+1)],[xn_(#+1),0]],color->[0,0.8,0],alpha->0.7);
connect([[xn_#,0],[xn_#,yn_#],[xn_(#+1),yn_(#+1)],[xn_(#+1),0]]);
sum=sum+dx*(yn_#+yn_(#+1))/2;
);
apply((-5..10),x,drawtext([x-0.1,-0.4],x,size->10));
drawtext([1,-1.5],"分割数を決めるスライダ",size->12);
drawtext([5,2],"分割数="+n,size->16);
drawtext([5,1.5],"台形の面積の和="+sum,size->16);
-----------------------------------------------------------------------------------------------
最初の関数の定義を変えればいろいろな関数について実験できます。
よくあるのが、
です。高校の教科書の例題にあります。x=tanθ とおく置換積分により、結果はπ/4になります。そこで、積分値を4倍して表示すれば円周率πの値が出るわけです。
シンプソンの公式
台形公式でもかなりよくなりますが、それでもまだ曲線との間に隙間がある。というわけで、上の辺を直線ではなく曲線にしてできるだけ隙間がないようにしましょう。そのときに、放物線を使うというのがシンプソンの公式です。
シンプソンの公式についてはWeb上を探せばいくらでも出てきますので、詳しくはそちらに譲るとして、ここではコンピュータ向けの（プログラム向けの）方法をとります。
シンプソンの公式の基本的な考え方はつぎの通りです。
まず、区間をn等分します。その一つの区間について、中点をとります。次の図は区間を2等分し、それぞれ中点をとったものです。薄い緑の線が中点のところの高さです。（区間を2n等分して、２つずつ使う、という考え方もあります）
この１つの小区間を[x0,x1,x2]とします。x1は区間[x0,x2]の中点です。
それぞれの関数値をy0,y1,y2としたとき、３点(x0,y0),(x1,y1),(x2,y2)を通る放物線で曲線を近似します。上の緑のラインがそれです。
この放物線の下方の面積を計算すると、
となります。（計算は定積分でおこないます。詳しくはWeb上で調べてください）
これを区間ごとに計算して、その総和をとればいいわけです。
さて、計算はこれで簡単に済みますが、それをシミュレーションとして表示するためには、近似した放物線を描かなければなりません。上の積分を行なうときには、放物線の方程式の各係数を求める必要はないのですが、描画するためには各係数を求めて、方程式を決定しなければなりません。高校数学では、点の座標を代入した３元１次連立方程式を解く問題なのですが、これを具体的な座標でなく、一般的な文字でおいた座標でやろうとすると結構面倒です。
これに関しては、いくつか方法があるようですが、ここではラグランジュ補間法というのを用いましょう。
ラグランジュ補間法を用いると、3点
を通る放物線の方程式は
となります。ラグランジュ補間法についてはWeb等で調べてください。ここで説明するにはスペースが足りません。
かくして、次のようなスクリプトを書けば、シンプソンの公式を使っての数値積分ができます。
-----------------------------------------------------------------------------------------------
f(x):=x^3-2*x^2+2;
A.y=0;
B.y=0; // 点A,Bがｘ軸から離れないようにする。
sx=A.x; //定義域左端
ex=B.x; //定義域右端
n=round((P.x-S.x)*2); // スライダから分点の数を得る
plot(f(#));
dx=(ex-sx)/n; //区間幅
xn=apply(0..n,sx+#*dx); // sxからexまでをn等分した分点のリスト
yn=apply(xn,f(#));
p(x,x0,y0,x1,y1,x2,y2):=(x-x1)*(x-x2)*y0/((x0-x1)*(x0-x2))+(x-x0)*(x-x2)*y1/((x1-x0)*(x1-x2))+(x-x0)*(x-x1)*y2/((x2-x0)*(x2-x1)); //3点を通る放物線
sum=0;
dx=dx/2;
repeat(n,i,
x0=xn_i;
x1=(xn_i+xn_(i+1))/2;
x2=xn_(i+1);
y0=yn_i;
y1=f(x1);
y2=yn_(i+1);
draw([x0,0],[x0,y0],color->[0,0.8,0],alpha->0.7);
draw([x1,0],[x1,y1],color->[0,0.8,0],alpha->0.7);
draw([x2,0],[x2,y2],color->[0,0.8,0],alpha->0.7);
plot(p(#,x0,y0,x1,y1,x2,y2),start->x0,stop->x2,color->[0,0.7,0]);
sum=sum+dx*(y0+4*y1+y2)/3;
);
apply((-5..10),x,drawtext([x-0.1,-0.4],x,size->10));
drawtext([1,-1.5],"分割数を決めるスライダ",size->12);
drawtext([5,2],"分割の数="+n,size->16);
drawtext([5,1.5],"近似した面積の和="+sum,size->16);
-----------------------------------------------------------------------------------------------
さきほどの台形公式での結果と比べてみましょう。この関数の区間[0,2]における定積分を計算すると
となりますから、確かにシンプソンの公式のほうが近い値になっています。2分割ですでにかなり近い値です。
また、
による円周率の計算もやってみましょう。さきほどの台形公式では12分割で3.1404でした。シンプソンの公式ではつぎのようになります。
2分割で3.1416と台形公式よりずっと近い値です。
なお、小数点以下の表示は4位までになっています。これはCindyScriptのdrawtextがそのように表示しているのですが、さらに細かい桁数で表示したい場合は、guess()関数を用いて
drawtext([5,1.5],"近似した面積の和="+guess(sum*4),size->16);
とします。
