1
大学数学基礎解説
文献あり

高校数学から始める数値計算(2):数値微分・数値積分

493
0

目次

本稿は「3. 数値微分」と「4. 数値積分」に当たります.なお,本稿はこの PDF版 を元に,Mathlogの仕様に合わせて一部の文言・体裁を変更したものです.内容は変わりません.

  1. はじめに
  2. 補間多項式
    1. ラグランジュ補間
    2. ニュートン補間
    3. 補間多項式による近似の誤差
  3. 数値微分
  4. 数値積分
    1. ニュートン・コーツの公式
    2. ガウス・ルジャンドル公式
  5. 補遺
    1. ルンゲ現象
    2. 補間多項式とテイラー多項式の関係

数値微分

関数を補間多項式で近似できるのなら,ある点における微分係数もまた,補間多項式の微分係数で近似できると考えられる.すなわち,実数x1,,xnはどの2つも相異なり,tx1,,xnに十分近いとき,点(x1,f(x1)),,(xn,f(xn))に関する補間多項式をp(x)とするとf(t)p(t)であると考えられる.

このような微分の近似を数値微分という.数値微分をするときはx1,,xnを等間隔に取ることが多い.以下にこの場合の例を示す.

数値微分の例1

n=2かつ(x1,x2)=(t,t+h)のとき
p(x)=f(t)+f[t,t+h](xt)
である.よって
f(t)dpdx(t)=f[t,t+h]=f(t+h)f(t)h
である.h0においては,この式の右辺は厳密にf(t)と等しくなる.

数値微分の例2

n=3かつ(x1,x2,x3)=(th,t,t+h)のとき
p(x)=f(th)+f[th,t](x(th))+f[th,t,t+h](x(th))(xt)
である.よって
dpdx(x)=f[th,t]+f[th,t,t+h](2x2t+h)f(t)f(t)f(th)h+h1(f(t)f(th))h1(f(t+h)f(t))(th)(t+h)h=f(t+h)f(th)2h
である.この式についても,h0においては
f(t+h)f(th)2h=f(t+h)f(t)2h+f(t)f(th)2hf(t)2+f(t)2=f(t)(h0)
が成立する.

数値微分の比較 数値微分の比較

図は,2つの直線
l1:y=f(t+h)f(t)h(xt)+f(t)l2:y=f(t+h)f(th)(t+h)(th)(x(th))+f(th)
f(x)=sin(x),t=π/3,h=0.2の場合について示したものである.図によれば,l2l1よりも傾きが1/2=cos(π/3)に近い.したがって,(sin(x))x=π/3における近似に関しては,例1よりも例2のほうが真値に近い値を与えている.より一般に,hを変えない場合,高次の補間多項式を利用した公式ほど近似値は真値に近づく傾向にある.

ニュートン・コーツの公式

関数f(x)の区間I=[a,b](a<b)における定積分を考える.実数x1,,xnIax1<<xnbを満たすとし,点(x1,f(x1)),,(xn,f(xn))に関する補間多項式をp(x)とおく.このとき,p(x)f(x)がよく近似できていれば
abf(x)dxabp(x)dx
となることが期待できる.

右辺について, 補間多項式の公式 を利用すると
abp(x)dx=ab(i=1nf(xi)kixxkxixk)dx=i=1nf(xi)abkixxkxixkdx
である.よって
wi=abkixxkxixkdx(i=1,,n)
とおくと,冒頭の近似式は
abf(x)dxi=1nwif(xi)
という形に書ける.

x1,,xnのことを標本点,実数w1,,wnのことを重みという.上の式から,重みw1,,wnは積分区間Iと標本点x1,,xnの取り方のみに依存し,関数f(x)によらないことが分かる.

特にここで,実数x1,,xnを区間[a,b]上で均等に取り
xi=a+(ba)i1n1(i=1,,n)
として定積分を近似する方法を閉じたニュートン・コーツの公式という(注:端点を含めず,xi=a+(ba)i/(n+1)とする方法を開いたニュートン・コーツの公式という.性質は閉じたニュートン・コーツの公式とほとんど同様である).

閉じたニュートン・コーツの公式の誤差は, 補間多項式による近似の誤差 の不等式を両辺ともに積分することで評価できる.実際
|abf(x)dxabp(x)dx|ab|f(x)p(x)|dx(ba)maxatb|f(t)p(t)|(ba)(maxaξb|f(n)(ξ)|n!(ba)n)=maxaξb|f(n)(ξ)|n!(ba)n+1
である(ただし,この評価は最良ではない.ニュートン・コーツの公式に関しては,より良い評価方法が知られている).

以下ではn=2,3について,w1,,wnを計算し誤差を見積もる.

台形則

n=2のとき
w1=abxx2x1x2dx=abxbabdx=ba2w2=abxx1x2x1dx=abxabadx=ba2
である.つまり,n=2のときの閉じたニュートン・コーツの公式は
abf(x)dxba2(f(a)+f(b))
と表される.この公式を台計則という.台計則の誤差は
|abf(x)dxabp(x)dx|maxaξb|f(2)(ξ)|2(ba)3
と評価できる.なお,台計則は2点(a,f(a)),(b,f(b))に関する補間多項式を利用しているので,1次以下の多項式関数に関しては( 補間多項式の一意性 により)厳密な値を与える.

なお,本稿では示さないが,もう少し厳しく
|abf(x)dxabp(x)dx|maxaξb|f(2)(ξ)|12(ba)3
と評価できることが知られている[3].

シンプソン則

n=3のとき
w1=ab(xx2)(xx3)(x1x2)(x1x3)dx=ab(x21(a+b))(xb)(a21(a+b))(ab)dx=ba6w2=ab(xx1)(xx3)(x2x1)(x2x3)dx=ab(xa)(xb)(21(a+b)a)(21(a+b)b)dx=2(ba)3w3=ab(xx1)(xx2)(x3x1)(x3x2)dx=ab(xa)(x21(a+b))(ba)(b21(a+b))dx=ba6
である.つまり,n=3のときの閉じたニュートン・コーツの公式は
abf(x)dxba6(f(a)+4f(a+b2)+f(b))
と表される.この公式をシンプソン則という.シンプソン則の誤差は
|abf(x)dxabp(x)dx|maxaξb|f(3)(ξ)|6(ba)4
と評価できる.

台計則と同様に考えれば,シンプソン則は2次以下の多項式関数について厳密な値を与えることが分かる.ところが,シンプソン則はさらに性質が良い.f(x)=x3のときを考えると,このとき
ba6(f(a)+4f(a+b2)+f(b))=ba6(a3+4(a+b2)3+b3)=b4a44=abf(x)dx
である.つまり,シンプソン則は3次多項式関数についても厳密な値を与える.

実は,シンプソン則の誤差は
|abf(x)dxabp(x)dx|maxaξb|f(4)(ξ)|2880(ba)5
と評価できることが知られている[3].この式からも,シンプソン則が3次以下の多項式関数について厳密な値を与えることが分かる.

台形則とシンプソン則の比較

台形則とシンプソン則 台形則とシンプソン則

図は,定積分31exdxの台形則とシンプソン則による近似で利用される補間多項式を示したものである.この場合,シンプソン則ではexの下側を通る部分と上側を通る部分がちょうど相殺され,より真値に近い値が得られると期待できる.実際,この真値は
31exdx=2.6685
であり,台形則とシンプソン則による近似値はそれぞれ
1(3)2(e3+e1)=5.53611(3)6(e3+4e1+e1)=2.8263
である.

閉じたニュートン・コーツの公式はたいてい,積分区間を複数の小区間に区切ってから,その各区間に対して適用される.この場合,台計則は複合台計則,シンプソン則は複合シンプソン則と呼ばれる(単に台形則,シンプソン則とも呼ばれる).例で見たように,閉じたニュートン・コーツの公式の誤差はbaに依存する.したがって,複合台形公式と複合シンプソン則の誤差は,baを小さく取るほど,すなわち,区間を多く分割するほど改良されると考えられる.

ガウス・ルジャンドル公式

ニュートン・コーツの公式では標本点を均等に取って定積分を近似した.この節では,標本点の取り方を工夫することで,ニュートン・コーツの公式と比べて少ない標本点で良い近似を得る方法を示す.

ルジャンドル多項式

正の整数nに対し
Ln(x)=12nn!dndxn[(x21)n]
で定義される多項式をn次のルジャンドル多項式という(注:正確には,この式におけるd/dxは形式微分である.本稿では普通の微分と特に区別せず扱う).

ルジャンドル多項式の例

L1(x)=122x=xL2(x)=18(12x24)=12(3x21)L3(x)=148(120x372x)=12(5x33x)L4(x)=18(63x570x3+15x)
である.

ルジャンドル多項式 ルジャンドル多項式

Ln(x)は区間I˚=(1,1)n個の相異なる根を持つ.

なお,多項式f(x)の根とはf(α)=0を満たす値αのことである.この定理を示す前に,次の補題を示す.

一般ライプニッツ則

関数f(x),g(x)n階微分可能とする.このとき,正の整数nに対して
dndxn(f(x)g(x))=k=0n(nk)f(k)(x)g(nk)(x)
である.ただし(nk)は二項係数である.

帰納法により示す.n=1のとき,積の微分法則により
ddx(f(x)g(x))=f(x)g(x)+f(x)g(x)
であり成り立つ.また,あるnでの成立を仮定すると
dn+1dxn+1(f(x)g(x))=ddx(dndxnf(x)g(x))=ddxk=0n(nk)f(k)(x)g(nk)(x)=k=0n(nk)(f(k)(x)g(nk+1)(x)+f(k+1)(x)g(nk)(x))=k=0n(nk)f(k)(x)g(n+1k)(x)+k=1n+1(nk1)f((k1)+1)(x)g(n(k1))(x)=f(x)g(n+1)(x)+(k=1n((nk)+(nk1))f(k)(x)g(n+1k)(x))+f(n+1)(x)g(x)=k=0n+1(n+1k)f(k)(x)g(n+1k)(x)((nk)+(nk1)=(n+1k))
となり,成立が確かめられる.

補題2により,定理1は以下のように証明できる.

xnk階微分したときの係数をAkとする.
f(x)=(x21)n=(x1)n(x+1)n
とおく.補題1により,l=1,,n1について
f(l)(x)=k=0l(lk)(dkdxk(x1)n)(dlkdxlk(x+1)n)=k=0l(lk)Ak(x1)nkAlk(x+1)n(lk)=(x1)nl(x+1)nlk=0l(lk)AkAlk(x1)lk(x+1)k
である.よって,f(1)(x),,f(n1)(x)x=1,1における値は0である.

このことを用い,l=1,,nについてf(l)(x)が区間I˚上にl個の相異なる根を持つことを示す.まず,l=1のときは平均値の定理により
f(ξ)=f(1)f(1)1(1)=001(1)=0
を満たすξI˚が存在する.

次に,あるl<nでの成立を仮定する.すなわち,l個の実数ξ1,,ξlI˚(ただしξ1<<ξlとする)が存在して,f(l)(ξ1)==f(l)(ξl)=0を満たすとする.このとき,ξ0=1,ξl+1=1とおくと,平均値の定理によりk=0,,lに対し
f(l+1)(ξk)=f(l)(ξk)f(l)(ξk+1)ξkξk+1=00ξkξk+1=0
を満たすξk(ξk,ξk+1)I˚が存在する.よって,f(l+1)(x)もまたl+1個の相異なる根をI˚上に持つ.

以上により,帰納的にf(n)(x)I˚上にn個の相異なる根を持つ.2nn!Ln(x)=f(n)(x)なので,Ln(x)もまたI˚上にn個の相異なる根を持つ.

k=1,,n1なら
11xkLn(x)dx=0
である.

定理1の証明と同様に
f(x)=(x21)n=(x1)n(x+1)n
とおく.このとき
11xkLn(x)dx=11xk12nn!f(n)(x)dx=12nn!11xkf(n)(x)dx
である.部分積分により
11xkf(n)(x)dx=[xkf(n1)(x)]1111kxk1f(n1)(x)dx
であるが,定理1の証明で見たように,l=1,,n1に対して
f(l)(1)=f(l)(1)=0
なので
11xkf(n)(x)dx=k11xk1f(n1)(x)dx
である.同様に部分積分を繰り返すことで
|11xkf(n)(x)dx|=|11k!x0f(nk)(x)dx|=k!|f(nk1)(1)f(nk1)(1)|=0
が得られる.よって
11xkLn(x)dx=12nn!11xkf(n)(x)dx=0
である.

以上により,次の定理がしたがう.

ガウス・ルジャンドル公式

I˚=(1,1)とし,x1,,xnI˚上にあるLn(x)の相異なるn個の根とする.このとき,f(x)2n1次以下の多項式関数であれば,点(x1,f(x1)),,(xn,f(xn))に関する補間多項式p(x)
11f(x)dx=11p(x)dx
を満たす.言い換えると,重みw1,,wn
wi=abkixxkxixkdx(i=1,,n)
で定めると
11f(x)dx=i=1nwif(xi)
が成り立つ.

Ln(x)の相異なるn個の根をx1,,xnとおく.関数f(x)に対し
G[f(x)]=i=1nwif(xi)(wi=abkixxkxixkdx)
とする.ニュートン・コーツの公式で見たように,G[f(x)]は点(x1,f(x1)),,(xn,f(xn))に関する補間多項式を区間[1,1]で定積分した値である.

G[]の定義により,Iで定義される任意の関数f(x),g(x)に関して
G[f(x)+g(x)]=G[f(x)]+G[g(x)]
が成立する.

f(x)2n1次以下の多項式関数とすると,f(x)は(多項式の意味で)Ln(x)により割れる.商をq(x),余りをr(x)とすると
f(x)=q(x)Ln(x)+r(x)
となる.したがって
G[f(x)]=G[q(x)Ln(x)]+G[r(x)]
であるが,x1,,xnLn(x)の根であるから
G[q(x)Ln(x)]=k=1nwkq(xk)Ln(xk)=k=1nwkq(xk)0=0
である.したがってG[f(x)]=G[r(x)]である.r(x)の次数はLn(x)の次数を超えないので,r(x)の次数は最大でもn1である.よって, 補間多項式の一意性 により点(x1,r(x1)),,(xn,r(xn))に関する補間多項式はr(x)に一致する.すなわち
G[f(x)]=G[r(x)]=11r(x)dx
である.

一方,q(x)の次数はn1以下なので,定理3により
11f(x)dx=11q(x)Ln(x)dx+11r(x)dx=11r(x)dx
である.よって
G[f(x)]=11r(x)dx=11f(x)dx
が成り立つ.

定理4により,関数f(x)2n1次以下の多項式関数でよく近似できるときは,標本点をルジャンドル多項式Ln(x)n個の根にすることで,より正確な近似値が得られると考えられる.この方法をガウス・ルジャンドル公式という.

ガウス・ルジャンドル公式とシンプソン則の比較

定積分31exdxをガウス・ルジャンドル公式によって近似する.

定積分31exdxは積分区間が[1,1]でない.そこで,置換積分により変数を
t=x+12,dx=2dt
に変え112e2t1dtを計算する.L3(x)の根x1,x2,x3(x1<x2<x3)は,方程式L3(x)=0を解くことで
12(5x33x)=0x1=35,x2=0,x3=35
と分かる.x1,x2,x3に対応する重みはそれぞれ
w1=59,w2=89,w3=59
である.したがって,ガウス・ルジャンドル公式によればf(x)=2e2x1とすると
112e2t1dt59f(35)+89f(0)+59f(35)=2.6651
である.これは例5で見たシンプソン則による近似値2.8263よりも真値2.6685に近い.

参考文献

[1]
堀之内 總一 and 酒井 幸吉 and 榎園 茂, Cによる数値計算法入門, 森北出版, 2015
[2]
金谷 健一, 数値で学ぶ計算と解析, 共立出版, 2010
[3]
菊地 文雄, 数値解析の原理: 現象の解明をめざして, 岩波数学叢書, 岩波書店, 2016
[4]
高橋 大輔, 数値計算, 理工系の基礎数学;8, 岩波書店, 2018
投稿日:202157
OptHub AI Competition

この記事を高評価した人

高評価したユーザはいません

この記事に送られたバッジ

バッジはありません。
バッチを贈って投稿者を応援しよう

バッチを贈ると投稿者に現金やAmazonのギフトカードが還元されます。

投稿者

数学、特に応用数学が好きです。Mathlogでは主に、数学とプログラミングを絡めたようなことを書けたらいいなと思っています。

コメント

他の人のコメント

コメントはありません。
読み込み中...
読み込み中
  1. 目次
  2. 数値微分
  3. ニュートン・コーツの公式
  4. ガウス・ルジャンドル公式
  5. 参考文献