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

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

参考までに、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になるなどの問題も発生する。


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

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



4次精度のルンゲ=クッタ法で問題を解いてみた [数値解析]

4次精度のルンゲ=クッタ法で問題を解いてみた


今日、11月24日に公開した数学の記事の第3問目は

  

になる。

この解は、問題3で求めたとおり

  

である。

(1)の連立常微分方程式は

  

の特殊なものであり、したがって、4次精度ののルンゲ=クッタ法で解くことができる。

連立常微分方程式を解く(4次精度の)ルンゲ=クッタ法のアルゴリズムは次のようなもの。

  

これを上から順次計算し、Step5まで計算しおえたら、xx+hだけ進め、Step1に戻って同じように計算すればよい。

ちなみに、上のアルゴリズムで使われている「=」は、数学の「等号」ではなく、右辺の計算結果を左辺に代入するという代入することをあらわしている。


4次精度のルンゲ=クッタ法はさすがに手計算でこれを行うのは大変なので、2次精度のルンゲ=クッタ法を使って、手計算で最初の値を求めてみることにする。


2次のルンゲ=クッタ法は次のようなアルゴリズム。

  


微分方程式(1)は

  

計算出発点のx=0のとき、y=1z=0

h=0.1とすると、

  


ちなみに、x=0.1のときの正しい値はy≒1.0050z≒0.10017だから、こちらも精度よく計算できていることがわかる。

そして、x=0.2のとき、上の計算で得られたy=1.0025z=0.1を使って同様の計算をする。

次に、4次精度のルンゲ=クッタ法を用いて、h=0.1としたときの計算結果を示す。




誤差は殆どなく、正確に計算できていることがわかると思う。


参考までに2次精度のルンゲ=クッタ法の計算結果を示す。




2次精度のルンゲ=クッタ法であってもこの程度の計算ならば、差は大きくない。

そこで、この両者の計算精度の差を明らかにするために、さらに複雑な次の連立常微分方程式を数値的に解くことにする

  

この微分方程式の解は、

  


4次精度の計算結果は次の通り。




厳密解とよく一致している。

しかし、2次精度のルンゲ=クッタ法では



と、xが大きくなるにつれて、厳密解との差が次第に大きくなってゆく。

h=0.1からh=0.01h1/10にしたとしても、精度はやはりh=0.1のときの4次精度のルンゲ=クッタ法に及ばない。



さらに1/10にしてh=0.001と、h=0.11/100したものが次のもの。


このとき初めて4次精度のルンゲ=クッタ法と2次精度のルンゲ=クッタ法は同程度の誤差になる。

精度が2次違うということはこういうことなんだケロよ。


4次精度のルンゲ=クッタのプログラムは以下の通り。

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

double f(double x, double y, double z) {
    return x*y*z; // 変えていいのはここ
}

double g(double x, double y, double z) {
    return x*y/z; //変えていいのはここ
}

// Runge_Kutta本体はいじらない方がいい
// いじらないでもいいように独立化させてある
// 独立化させているので、この部分は他のプログラムにそのまま流用可能!!
void Runge_Kutta(double x, double y, double z, double h ,double n) {
    double dely1, dely2, dely3, dely4;
    double delz1, delz2, delz3, delz4;
    int i;

    printf("%f %f %f\n",x,y,z); // 結果の出力(初期値)
    for (i = 1; i <= n; i++) {
        dely1 = h*f(x,y,z);
        delz1 = h*g(x,y,z);
        dely2 = h*f(x+h/2, y+dely1/2, z+delz1/2);
        delz2 = h*g(x+h/2, y+dely1/2, z+delz1/2);
        dely3 = h*f(x+h/2, y+dely2/2, z+delz2/2);
        delz3 = h*g(x+h/2, y+dely2/2, z+delz2/2);
        dely4 = h*f(x+h, y+dely3, z+delz3);
        delz4 = h*g(x+h, y+dely3, z+delz3);
        x=x+h;
        y=y+(dely1+2*dely2+2*dely3+dely4)/6.;
        z=z+(delz1+2*delz2+2*delz3+delz4)/6.;
        printf("%f %f %f\n",x,y,z); // 結果の出力
    }

}

main() {
    int n;
    double x, y, z;
    double h;

    h=0.1; // xの増分の設定 ココは変えてもよい
    n=10;  // 何回計算するかの設定 最終的にn×h進む
    x=1; y=1./3.; z=1; // x,y,zの初期条件の設定
    Runge_Kutta(x,y,z,h,n); // Runge-Kutta法を呼び出す
   
    /* dy/dx=f(x,y,z), dz/dx=g(x,y,z)の初期値問題を
       4次精度のルンゲ=クッタ法で解くプログラム */
    // 注意 このプログラムで解く関数y,zはx=√7で不連続、ココでゼロ割発生!!
}


上のプログラムは、2つ目の微分方程式を解いている。
1つ目の微分方程式を解きたい場合、

double f(double x, double y, double z) {
    return z; // 変えていいのはここ
}

double g(double x, double y, double z) {
    return y; //変えていいのはここ
}

とし、mainの
    x=1; y=1./3.; z=1; // x,y,zの初期条件の設定

 x=0; y=1; z=0;
にすればよい。

陽解法による拡散方程式の解法の安定性 [数値解析]

陽解法による拡散方程式の解法の安定性


次の拡散方程式がある。

  

ここで、uは時刻tと位置xの関数、すなわち

  

とする。

数値計算の手法にならって

  

と書くことにする。

陽解法は(1)を次の差分方程式に置き換えて、次の差分方程式を解くことによって(1)の近似解を求める。

  

とおくと、(2)は

  

数値解と厳密解の誤差は

  誤差=数値解−厳密解

である。

数値解、厳密解ともに(3)を満たすので、誤差も(3)を満たす。

つまり、

  

誤差が次のように表せるとすると、

  

(4)より

  

オイラーの公式より

  

だから、

  

さらに、半角公式

  

より

  

で、

  

だから、 (5)式の右辺は誤差の倍率をあらわすと考えられる。

したがって、

  

のときに、は収束する。

で、

  

とおくと、

  

を満たすK

  

したがって、

  

これをすべてのkについて満たさなければならないから、

  


上の議論は正確なものではないので注意。

このようにして、陽解法の

  

という制限、条件が出てくるという話です。

さてさて、この話は本当かということで、この検証のためのスプレッドシートを作り、実験してみた。

解く偏微分方程式は、

  

これをΔt=0.25Δx=0.5で解いた結果はこれ。
x=1とx=2の計算結果を示してある。




厳密解、そして、陰解法で解いた結果と良好な一致を見せており、このスプレッドシートの正当性を示している。


r=0.5
の場合





r=0.6
の場合。



時間の経過とともに誤差が次々と伝播、増大し、数値解が激しく振動してしまう。




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