So-net無料ブログ作成
検索選択
数値解析 ブログトップ
前の10件 | -

オイラー法とシンプレクティック法を用いた1次元調和振動子の解法 [数値解析]

オイラー法とシンプレクティック法を用いた1次元調和振動子の解法

 

1次元調和振動子のハミルトニアン(力学的なエネルギー)は

  

である。

ここで、

  

であり、qはバネの自然長からの変位、pは運動量、(1)式のp²/2は運動エネルギー、q²/2はバネの位置エネルギー(ポテンシャルエネルギー)に相当する。

 

(2)式は、

  

と書き換えることが可能で、(2)式とバネの単振動に関するニュートンの運動方程式と同一のものである。

時刻t=0における初期条件を(q,p)=(1,0)とすれば、(2)、(3)式の解は

  

である。

 

1次精度のオイラー法を用いて解くには、

  

つまり、

  

と近似し、これを逐次繰り返して計算すればよい。

 

たとえば、t=0のとき、(q,p)=(1,0)であれば、Δt=0.1とすれば、

  

この計算結果をもとに、

  

といった具合に計算できる。

 

この計算は非常にシンプルな繰り返し計算だから、プログラムを作らずとも、表計算ソフトで簡単に計算が可能で、以下にt=0における初期条件を(q,p)=(1,0)としたものの計算結果を示す。




Euler-hyou-graph.png
 

このように非常に粗い近似に基づく計算でありながら、比較的によく計算できていることがわかるだろう。時間の経過とともにqpの厳密解と1次精度のオイラー法を用いた数値解との食い違いが大きくなるのは、近似計算を繰り返すたびに誤差が伝播し、蓄積するため。このため、オイラー法の誤差は雪だるま式に増加し、計算はやがて発散する。

次に、オイラー法を用い、Δt=0.1とし、t=10まで数値計算結果を(q,p)のグラフとして表したものを以下に示す。


Euler1st.png
 

解いた問題は、力学的エネルギー、つまり、ハミルトニアンHが一定であり、

  

(q,p)は単位円上に存在しなければならないのだが、時間の経過とともに軌道q²+p²=1からの逸脱が大きくなってゆく。

つまり、オイラー法を用いた近似計算法はエネルギー保存則を満たす解法ではなないということだ。

 

オイラー法のこの欠点を改良したものにシンプレクティック(Sympretic)法と呼ばれるものがある。

オイラー法で用いる計算式(4)を次のように改良したもの。

 

  

 

この僅かな改良によってエネルギー保存則からの逸脱がかなり改善され、オイラー法の宿命的な欠点である時間の経過に伴う誤差の伝播、蓄積が解消される。

 

以下に、表計算ソフトを用い、シンプレティック法を用いて解いた計算結果を示す。




Symplectic-hyo-graph.png 

 

比較参照のために、シンプレクティク法とオイラー法を用いて解いたqと厳密解をのグラフを以下に示す。


Chowa-Euler-Symplectic-graph-001.png
 

さらに、シンプレティック法による数値解を(q,p)のグラフで表したものを示す。


Symplectic.png
 

オイラー法とは異なり、シンプレクティック法を用いた数値解はかなり良くq²+p²=1の円周上に存在していることがわかるだろう。

シンプレクティク法とオイラー法を用いて解いたq,pを元に計算したハミルトニアンHと厳密な値1/2との相対誤差

  

の時刻による推移をグラフに示す。


gosa.png
 

このグラフが示すように、1次精度のシンプレクティク法は完全にエネルギー保存則を満たすわけでないが、振動をしつつも常にある誤差の範囲内に収まり、オイラー法のように指数関数的に誤差が増大することはない。

参考までに2次精度をのシンプレクティク方による計算結果も示してある。

 

数値計算ではよくあることなのだけれど、(4)式を少し変え(5)式にし計算するだけで計算結果が劇的に変化するのだから、数値計算は奥が深いよね。

そう思いませんか?

そして、少しでもこうした数値計算に興味を持ってもらえれば幸せです。

 

なお、これらの議論をより正確にするためには、ニュートン力学の進化形である解析力学などの知識が必要になるので、感覚的に理解できれば十分ではないでしょうか。

 

参考までに2次精度のシンプレクティック法で解いた計算結果は以下に示す。


Symplectic2nd.png


ネムネコ、シンプレクティック法を使って微分方程式の数値積分をする!! [数値解析]

かつてのネムネコ・ファミリー――某Q&Aサイト上のみでネムネコと深い交流のあったヒトたち――であるddt³さんから、シンプレクティック(Symplectic)法と呼ばれる数値計算法、微分方程式の数値積分法を教えてもらったので、バタバタとプログラムを作って計算してみた。

 

数値計算につかったモデルは、調和振動子

  

で、微分方程式は

  

で、t=0のときの初期条件を(q,p)=(1,0)として、1次精度のオイラー法と1次精度のシンプレクティック法を使って、この微分方程式を解いてみた。

進めるタイムステップΔt=0.1として数値的に解いた。

 

(1)式のHはハミルトニアンと呼ばれるもので、この場合、運動エネルギーp²/2とバネの位置エネルギーq²/2の和になる。力学的エネルギーが一定だから、(2)を数値的に解いた(q,p)をプロットすると、本来、数値解は

の円周上に存在しなければならない。

しかし、1次精度のオイラー法を用いて解いた数値解は時間の経過とともにp²+q²=1の円から離れていくんだケロよ。つまり、エネルギー保存則が満たされていない!!

軌道からどんどん離れてゆくから、時間の経過とともにエネルギーが増えていく!!


Euler1st.png
1次精度のオイラー法を用いた数値計算結果


これに対し、シンプレクティック法は、多少振動するけれど、力学的なエネルギーが保存されるんだケロ。


Symplectic.png
1次精度のシンプレクティック法を用いた計算結果


数値積分法は、

1次精度のオイラー法だと

  

1次精度のシンプレクティック法だと

  

だけの差なのだけれど、この僅かな差でこれほど計算結果が異なるというのは正直驚きだにゃ。



ネムネコは、日々進化しているにゃ。




ニュートン法のプログラム [数値解析]

参考までに、C言語で書いたニュートン法のサンプルプログラムをを以下に示す。

#include <stdio.h>
#include <math.h>

double f(double x) {
    return tan(x)-1.;
}

double df(double x) {
    return 1/(cos(x)*cos(x));
}

double Newton(double x, double eps, int n) {
    double x0;
    int i;

    for (i = 0; i < n; i++) {
        x0=x;
        x=x-f(x)/df(x);
        if (fabs(x-x0) < eps ) return x;
    }
   
    return -1;
}

// f'(x)の計算で中心差分を用いたニュートン法
double Newton2(double x, double eps, int n) {
    double x0, h = 0.001;
    double df;
    int i;

    for (i = 0; i < n; i++) {
        x0=x;
        df=(f(x+h)-f(x-h))/(2.*h); // f'(x)の計算に中心差分を使う
        x=x-f(x)/df;
        if (fabs(x-x0) < eps ) return x;
    }
   
    return -1;
}   

main() {
    double x=0.0, eps =1.e-6;
    int n = 10;

//    f(x)=tan x  1=0を解く
    printf("x=%f\n",Newton(x,eps,n)); // ニュートン法
    x=0.0;
    printf("x=%f\n",Newton2(x,eps,n)); // f'(x)に中心差分を用いたニュートン法
}

収束判定として、|x−x₀|<10⁻⁶、または、計算の反復回数が10を越えたときを用いている。

解く方程式はtan x = 1で、0≦x<π/2とすると、この解は
 x=π/4≒0.785398163

反復回数6回で|x−x₀|<10⁻⁶を満たし、
 x=0.78598
という収束解に到達する。

上のプログラムでは、f'(x)の値を中心差分を用いて計算するニュートン法のプログラムもあわせて紹介している。
中心差分を用いると、微分係数f'(x)は

  

と近似できることを利用し、この値を実際の微分係数f'(x)の値のかわりに使用している。
つまり、中心差分を利用したニュートン法の場合、導関数f'(x)を与える必要がなく、f(x)だけをプログラムに記述すればよい。

このプログラムでは、あえてh=0.001と、hに粗い値を与えているが、微分係数と中心差分を用いて微分係数との誤差はおおよそh²=0.001²=10⁻⁶のオーダーなので、このように粗いhの値でも
 x=0.78598
という解が得られており、実用上、問題がないことがわかると思う。


ニュートン法 [数値解析]

ニュートン法

Newton-method-01.pngf(x)を微分可能な関数とする。

方程式

  

の解をx=αとする。

右の図のように適当な点x₁を選び、y=f(x)の点(x₁,y₁)における接線の方程式は

  

で、この接線とx軸との交点のx座標x₂は、上の式にy=0を代入することによって、

  

となる。

そして、同様ににおける接線の接線を引き、この接線とx軸との交点のx座標x₃を求めると、

  

となり、この操作を繰り返せば繰り返すほど、こうして求められたf(x)=0の近似値である

f(x)=0の解、x=αに近づいてゆくことが予想される。


これがニュートン法である。

漸化式の形で書けば、ニュートン法は次のようになる。

  


Newton-method-tab-01.pngx=√3の近似値をニュートン法を用いて求めることにする。

x=√3の両辺を2乗すると

  

これを

  

とおき、f(x)=x²−3=0x>0)とすれば、これはx=√3と同値。

f'(x)=2xだから、

  

計算開始のx₀=1として、表計算ソフトを使って計算したものは次の通り。

4、5回計算するだけで、x=√3≒1.732050808という近似値に到達している。


ニュートン法は前回の2分法よりも速く、しかも急速に収束することがわかると思う。


ただし、ニュートン法は、次の例のように、収束しないことがある。

  

計算の初期値としてx₀=1、または、x₀=2を取ると、

接線の方程式が、それぞれ、

  

となり、x=1x=2を交互に永遠に行き来する。

Newtonhou-02.png

この他にも、f'(x)=0になる点に差し掛かったとき、ゼロ割が発生するなど、危険な一面も有している。

こういうことは極まれにしか起きないけれど、運悪くこのような事態に遭遇することがある。

2分法と比較すると、ニュートン法は収束の速度は速いけれど、安定性に欠ける。


前進差分と中心差分 [数値解析]

前進差分と中心差分

2次関数f(x)=ax²+bx+ca≠0)があるとする。

この導関数f'(x)

  

だから、x=αにおける微分係数f'(α)

  

である。

また、このとき、x=αにおける中心差分は

  

となる。

したがって、中心差分は、f(x)が1次関数、2次関数の場合、微分係数を正確に計算する。


f(x)=sinxとし、[−π,π]を20分割し、前身差分、中心差分によってf'(x)の近似値を求めたものを以下の図に示す。


  

と粗い計算であるが、中心差分の方が前進差分よりも精度よく計算出来るだけでなく、f'(x)=cosxを再現していることがわかると思う。


微分と差分 [数値解析]

微分と差分


f(x)を何度でも微分可能な関数とする。

このとき、f(a+h)は、以下のようにテーラー展開が可能である。

  

同様に

  

(1)式から

  

h1次以上の項を無視すると、x=aにおける微分係数f'(a)は次のように近似可能である。

  

同様に、(2)式から

  

が得られる。

(3)、(4)式とも1次以降の項を落としているので、誤差は1次オーダ、すなわち、O(h)である。

また、(1)から(2)を引くと

  

以上の項を落とすと、

  

という近似式が得られる。

(5)式は、以降の項を落としているので、誤差はhの2次オーダ、O(h²)

x軸上に等間隔hでならぶの点の集まりに正の方向に向かって整数の番号をつけると、という点の列、点の列が得られる。

とすると、(3)、(4)、(5)式は次のように書き換えることができる。

  

このように書くのは面倒なので、と略記することにすると、(3)、(4)、(5)の近似式は

  

となる。

このように、無限小の微分を有限な差の形で近似する方法を差分法という。

(6)を前進差分、(7)を後退差分、(8)を中心差分という。

2次導関数f''(x)の近似式は、次のように求めることができる。

(1)と(2)を足すと、

  

以上の項を落とすと、

  

つまり、

  

となる。

なお、(6)、(7)式は不等間隔で点がならんでいる場合でも、(6)式では、(7)式はとおくことによって成り立つが、(8)、(9)は不等間隔の場合、成り立たないので注意が必要。

不等間隔の場合、テーラー展開にさかのぼって、差分による微分の近似式を求めないといけない。

x=0における微分係数を(6)、(7)、(8)を使って求めてみる。

h=0.1である。

  

前進差分、後退差分の誤差は約0.05、中心差分は約0.0017だから、この場合、中心差分が最も精度よく計算できていることがわかるだろう。

sabun-gosa-01.pngx=2のときの微分係数の誤差とxの増分hとの関係を右図に示す。

前進、後退差分の勾配が1、中心差分の勾配が2であることがこの図から分かると思う。

このことは、前進・後退差分の誤差のオーダーが1次であり、中心差分の誤差のオーダーが2次であることを指し示しており、理論通りというわけ。
つまり、h1/10になったとき、前進、後退差分の精度が10倍良くなるのになるのに対して、中心差分は精度は100倍向上する。



数値積分 台形公式、中点公式とシンプソンの公式の導出と誤差 [数値解析]

数値積分 台形公式、中点公式とシンプソンの公式の導出と誤差


次の定積分を考える。

  

これは、x=t+aと変数の変換を行えば、dx=dtでかつ、x=aにはt=0x=a+hにはt=hが対応するので、

  

になる。

関数f(x)が何度でも微分可能、つまり、級であるとき、f(a+t)は次のようにテーラー展開することが可能。

  sekibun-gosa-01.png

あるいは、

  

(2)を(1)に代入すると

  

総和記号Σの中は
  sekibun-gosa-03.png

したがって、

  

ここで、

  

と前進差分を使うと、
  

a+h=bと置けば

  sekibun-gosa-05.png

したがって、

  

と、それを台形で近似した

  

との誤差は程度ということになる。


そして、この結果を用いて、

  

という台形公式の誤差の限界公式を導くことができる。

より厳密な議論は、たとえば、ねこ騙し数学の次などを見て欲しい。


台形公式の精度を求める問題


積分区間[a,a+h]の中点をcとし、x=t+cという変換をすると、dx=dtで、x=aにはt=−h/2x=a+hにはt=h/2が対応するので、

  

f(c+t)をテーラー展開すると

  

(4)に(5)を代入すると、nが奇数のとき

  

nが偶数のとき

  

になるので、
  

になる。

sekibun-gosa-fig-01.pnga+h=bとおくと、

  

だから、

  

そして、

  

と近似する方法を中点公式と呼ばれる。


上の議論から、中点公式の誤差は

  

程度で、誤差のオーダーはである。

[a,b]f(x)>0のとき,

  

は長方形DEFGの面積。

また、

  sekibun-gosa-07.png

として、(6)式に代入すると、

  

これから、シンプソンの公式の誤差が

  

程度で、シンプソンの公式の誤差がh⁵のオーダーであることが分かる。

シンプソンの公式は

  sekibun-gosa-09.png

で、(8)式と

  

とでは、右辺の分母が36と異なっているが、これは(8)式ではb−a=2hとしているのに対して、(9)式ではb−a=hとしていることに由来する。

(9)式の形で(8)式を書きなおすとき、h2hに変えれば良いので、次のようになる。

  

2次精度のルンゲ=クッタ法 [数値解析]

2次精度のルンゲ=クッタ法


αβを適当に選ぶと、所定の誤差の範囲で、における定積分を

  

と近似できるものとする。

微分方程式は

  

よって、

  


ところで、を点でテーラー展開し、で打ち切ると

  

また、

  

だから

  

よって、

  


一方、同じくテーラー展開し、1次で打ち切ると

  

(2)と(3)を(1)に代入すると、

  

だから、

  


以上のことまとめると、

  


これが2次精度のルンゲ=クッタ法と呼ばれるものであり、その導出の流れるになる。


これをコンピュータや手計算で近似計算する場合、次のようにすればよい。


2次精度のルンゲ=クッタ法

  


修正オイラー法は

  


2nd_R-K_VS_IE.png2次精度のルンゲ=クッタ法と修正オイラー法を用いて次の微分方程式を数値的に解いた結果は次の通り

  
この結果から、2次精度のルンゲ=クッタ法と修正オイラー法の精度が同程度であることが理解できると思う。




修正オイラー法 [数値解析]

修正オイラー法


1階の常微分方程式の初期値問題の一般形は

  

だが、より簡単な次のものを考える。

  


s-graph-01.png(2)の解は

  

g(x)>0の場合、(3)の右辺

  

g(x)x軸と2直線で囲まれた部分の面積Sに等しい。


オイラー法の場合、(4)の積分は

  

で近似される。

右図で正しい面積はピンク色の部分、(5)式で与えられる面積は斜線部で示されている。

そして、このとき、(3)は

  


s-graph-02.png一方、式(4)のこの定積分を台形公式で近似すると

  

(右図参照)

この図から明らかなように、(5)より(7)の方がより面積を精度よく近似しており、(3)を

  

と近似したものの方が精度よく計算できるものと考えられる。

(1)と(2)はまったく同じ問題ではないけれど、yxの関数であり、したがって、f(x,y)xの関数であり、f(x,y)=g(x)と考えれば、(8)の結果を使って、

  

が得られる

この近似計算法が修正オイラー法と呼ばれるものである。


⑨の右辺にが出てくるが、これは左辺のと、本来、同一のもので未知数。したがって、⑨を使ってを求めるためには、何らかの方法で右辺のを推測する必要がある。
そこで、

   

と近似することによって、⑨の右辺を計算することができる。


修正オイラー法をもちいて微分方程式の初期値問題を近似的に計算するには次のようにすればよい。

  


微分方程式

  

を、分割幅h=0.1として、x=0からx=3まで、オイラー法と修正オイラー法を用いて数値的に解いたものが下図である。

s-graph-03.png

赤の+はオイラー法、青の*が修正オイラー法を使って求めた数値解、そして、実線がこの微分方程式の解、すなわち、

  

である。

図を見れば分かるように、修正オイラー法の方が厳密解とよく一致している。


seido-Euler-graph.png計算する区間[a,b]の分割数nと計算精度の関係を右の図に示す。

精度の比較のために使用した微分方程式は

  

a=0b=2とし、計算する領域を[0,2]としている。

このグラフは、縦軸に厳密解との誤差、横軸に分割数をとったもので、縦軸、横軸ともに対数目盛りを使っている。

このグラフを見ると、オイラー法では分割数が10倍になると誤差がおよそ1/10に減少し、修正オイラー法では分割数が10倍になると誤差が1/10²になっていることがわかると思う。
つまり、オイラー法の計算精度が分割数nの逆数(分割幅h)に比例しているのに対し、修正オイラー法は分割数の逆数(分割幅h)の2乗に逆比例するのであった。

参考までに、この計算に使用したプログラムを以下に示す。

#include <stdio.h>
#include <math.h>
#define N 100000

double f(double x, double y) { // dy/dx=f(x,y)=x-y
    return x-y;
}

double exact(double x) { // dy/dx=x-yの厳密解
    return 2.*exp(-x)+x-1;
}


// オイラー法
void Euler(double x, double *y, double h, int n) {
    int i;
   
    for (i=0; i <n; i++) {
        y[i+1] = y[i]+f(x,y[i])*h;
        x=x+h;
    }
}

// 修正オイラー法
void IEuler(double x, double *y, double h, int n) {
    double dy1, dy2;
    int i;
   
    for (i=0; i <n; i++) {
        dy1 = h*f(x,y[i]);
        dy2 = h*f(x+h,y[i]+dy1);
        y[i+1] = y[i]+0.5*(dy1+dy2);
        x=x+h;
    }
}


main() {
    double y[N+1];
    int n[] = {10, 50, 100, 500, 1000, 5000,10000};
    int i;
    double x0, y0,h;
    double a,b;
    double x;
    double err1, err2;

    a=0; b = 2.;
    x0=a; y0=1.;
    for (i=0; i < 7; i++) {
        h=(b-a)/n[i];
        x = x0; y[0] =y0;
        Euler(x0,y,h,n[i]);
        err1=fabs(y[n[i]]-exact(b));
        IEuler(x0,y,h,n[i]);
        err2=fabs(y[n[i]]-exact(b));
        printf("%d %f %12.8f %12.8f1 \n",n[i],h,err1, err2);
    }
}

多項式による補間 [数値解析]

多項式による補間


§1 補間


2つの点の間に位置する点、すなわち、におけるf(x)の値を求めることを補間するという。2つの点の外の点におけるf(x)の値を求めることを補外するという。


もちろん、例えば、f(x)=x²+3x+2のように関数が式の形で与えられ、しかも、この計算が簡単にできる場合、補間や補外をする必要はなく、f(x)=x²+3x+2を用いて計算すればよい。しかし、f(x)=sinxならば、電卓などを使わない場合、値が既知の点をもとにテーラー展開などをして値を求めざるを得ない。

また、複数の点、たとえば、の値は与えられているけれど、y=f(x)の関数が与えられていないことも多く、何らかの手段で既知でない点の値を求めざるを得ない。

このような時に、補間や補外といった手法が使われる。

補間や補外には、テーラー級数がよく用いられる。すなわち、

  

である。

と書くのは面倒なので、と書くことにすると、

  

これを2次の項、すなわち、で打ち切ると、

  

に前進差分を

  

を用いると、

  

になる。

ここで、

  

とおくと、(2)は

  

となる。

この(3)を1次補間公式という。


なにか難しいことを書いているように思うかもしれないが、(3)は、における値f(x)を2点を通る直線(線分)で近似したものに過ぎない。


(1)を3次の項で打ち切るときには、

  

として、これを(1)に代入すると、
  

という2次補間公式が得られる。


例1 3点(0,1)(1,−1),(2,3)があるとする。

  

とすると、
  

この結果を(1)に代入すると、

  

という補間曲線が得られる。

lagrange-graph-01.pngそして、これが2次補間によって得られる曲線である。


もちろん、

  

として、f(0)=1f(1)=−1f(2)=3としてえられる次の連立方程式を解いてもよい。

  

この連立方程式の解は(a,b,c)=(3,−5,1)となり、同じ結果が得られる。

 


§2 ラグランジュ補間


次のn次の多項式を考える。

  

この多項式は以外のすべての点

  

において0である。

そこで、

  

とおくと、多項式において1となり、それ以外の0になる。

つまり、

  

だから、

  

とおくと、このn次の多項式はn+1個の点

  

を通る多項式になる。


これがラグランジュ補間と呼ばれる方法である。


例2 3点(0,1)(1,−1),(2,3)があるとする。

  lagrange-04.png

だから、

  



lagrange-tab-01.png
 上の表の値をもとに、ラグランジュ多項式によってf(6)を補間せよ。

【解】

  lagrange-06.png

よって、

  

(解答終了)

これが人間向きの計算法かと言えば大いに疑問であるが、単純な繰り返し計算を得意とするコンピュータにとっては都合のいい計算法である。


以下にこのサンプルプログラム(C言語)を示す。


#include <stdio.h>
#include <math.h>

double f(double x) {
    return sin(x); // f(x)=sin x   [−π,π]
}

// ラグランジュ補間
double lagrange(double *x, double *y, int n, double t) {
    double p,s;
    int i,j;
   
    s=0;
    for (i = 0; i <= n; i++) {
        p = 1.;
        for (j=0;j<=n; j++) {
            if (i != j) {
                p = p*(t-x[j])/(x[i]-x[j]);
            }
        }
        s = s+y[i]*p;
    }
    return s;
}

main() {
    double x[10]={0.0};
    double y[10]={0.0};
    double t,h;
    int i;
    int n = 9; //n次式
    double a, b;

    b = 4.*atan(1); // 区間[a,b]のbの値πをセット
    a = -4.*atan(1); // 区間[a,b]のaの値−πをセット
    h = (b-a)/n; // [a,b]をn等分した幅
    for (i = 0; i <= n; i++) { // i=0〜nの点の値をセット
        t = a+i*h;
        x[i]= t;
        y[i] = f(t);
    }

// 計算結果の出力
    printf("    x      lagrange   sin x \n");
    for (i=0; i <=100; i++) {
        t=a+i*(b-a)/100;
        printf("%f %f %f\n",t, lagrange(x,y,n,t), f(t));
    }
}


f(x)=sinx(−π≦x≦π)を5次のラグランジュ補間で近似計算した値を以下に示す。この場合、よく近似できていることがわかる。

lagrange-graph-02.png




誤解がないように言っておきますが、与えられた点をすべて通る多項式を求めただけで、ラグランジュ補間などの(等間隔の)多項式補間が元の曲線を忠実に反映しているというわけではない。


lagrange-graph-03.png例えば、

  

という曲線があるとする。

これを[−1,1]で等間隔に5分割(点の数は6)、9分割(点の数は10)し、その値をもとに5次補間、9次補間すると右のような曲線が得られる。


点が多く、高次で補間するほうがよりもとの曲線を反映するように考えられるが、この関数の場合、高次にすればするほど、端点であるx=±1の近傍でもとの曲線との乖離が大きくなる現象、ルンゲ現象が発生する。


直線で結んだほうがこの曲線の挙動を正確に反映している。


ここで言いたいのは、データがn+1個あるからn次の多項式で近似するということは非常に危険な行為だということです。


n次のラグランジュ多項式で結べば確かにのすべての点を通る曲線を得られるけれど、時に、このルンゲ現象のようなことが発生することがある。


また、この曲線の極値はx=0のところにのみあるにもかかわらず、5次補間では極値をとる点がが3、9次補間では極値をとる点が7になるなどの問題も発生する。


だから、物理や化学などの実験で得られたデータなどのグラフは、絶対に高次のラグランジュ補間で結んではいけない。

実験データを高次のラグランジュ補間で結んだグラフを書き、それをレポートなどにのせて提出すると、必ず、実験担当の先生や指導教官からこっぴどく叱られる。「このデータを曲線で結んだ理論的根拠は何だ。説明しろ!!」と問い詰められるのが落ちである。



前の10件 | - 数値解析 ブログトップ