ホイン（Heun）法

ホイン（Heun）法

 

§１　オイラー法

 

微分方程式

の解は

である。

（２）式の右辺の積分を

で近似すると、（２）式は

と近似することにより、x₀を計算の起点とし、x₁、x₂、・・・と築時的に計算をすることができる。

とあらわすことにすると、

になる。

とくに、点が等間隔hに並んでいるとき、

である。

この方法をオイラー法という。

 

問題　オイラー法を用いて、次の微分方程式のx=0.1、0.2、0.3におけるyの値を求めよ。

【解】

x₀=0、y₀=y(x₀)=1、h=0.1とおくと、

x₁=0.1、y₁=1.05だから

x₂=0.2、y₂=1.163だから

したがって、

x=0.1のとき、y=1.05

x=0.2のとき、y=1.1106

x=0.3のとき、y=1.1846

（解答終）

 

表計算ソフトを用いた計算結果を以下に示す。

 

 

 

オイラー法だと、計算を進めると、誤差が蓄積して、その結果、厳密解

との食い違いが大きくなっていくことが理解できると思う。

 

 

§２　ホイン法

まず、オイラー法

を用いてを計算し、次に台形公式

と修正する。

以下同様に、所定の精度が得られるまで

と反復計算をする方法をホイン法という。

 

ホイン法を用いて、問題の微分方程式をh=0.1と同一の条件で解いた結果は次の通り。

 

 

 

 

この計算に用いたプログラムは以下の通り。

 

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

double f(double x, double y) {
    return 0.5*(1+x)*y*y;
}

void Heun(double *x , double *y , double h , int n) {
    double eps = 1.0e-6, err;
    double yn = 0.;
    int i,j, limit = 100;

    for(i=0; i< n; i++) {
        x[i+1]=x[i]+h;
        // Euler法で予備計算
        y[i+1] = y[i] + h*f(x[i],y[i]);
        for (j=1; j <= limit; j++) { // 所定の精度が得られるまで反復
            // 台形公式で修正
            yn = y[i]+0.5*(f(x[i],y[i])+f(x[i+1],y[i+1]))*h;
            err =fabs(y[i+1]-yn);
            if (err <= eps) break;
            // y[i+1]の更新
            y[i+1]=yn;
        }

    }
}

double Exact(double x) {
    return 4./(4 - 2*x -x*x);
}

main() {
    double x[11], y[11];
    double h;
    int i, n;

    h=0.1;
    n=10;

    for (i = 0; i<=n; i++) {
        x[i] = 0; y[i]=0;
    }
   
    x[0]=0; y[0]=1.; // xとyの初期値
   
    Heun(x,y,h, n); // Heun法を呼び出す

    printf("   x      y(Heun)  y(exact) error(％) \n");
    for (i=0; i <= n; i++) {
        printf("%f %f %f %f\n", x[i],y[i],Exact(x[i]),fabs(1.0-y[i]/Exact(x[i]))*100);
    }
}

 

オイラー法より、Heun（予測修正子法ともいう）のほうが精度よく計算できていることがわかるだろう。

タグ：数値解析

数値計算の誤差：打ち切り誤差と丸め誤差

数値計算の誤差：打ち切り誤差と丸め誤差

 

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

このとき、f(a+h)はx=aにおいて次のようにテーラー展開が可能である。

有限項で打ち切れば、

は剰余項であり、有限項で打ち切ったときの誤差に相当する。

 

最も簡単なa=0のときを考えてみる。

かりに

という関数があるとする、

だから、これをx=0でテーラー展開（マクローリン展開）すると、

と一致する。

 

有限な多項式で表される関数ならば、上の例のように第n次以下の導関数を求めることによってその係数を求めることができる。しかし、テーラー展開（マクローリン展開）が指数関数のように、

と無限級数になる場合、これを無限に計算することはできないので、適当な次数nで打ち切って計算をしなければならない。

たとえば、３次の以上の項を落とせば、

そして、その誤差、打ち切り誤差O(x³)は

と評価されるであろう。

 

この他に、有限桁の数しか扱えないコンピュータの計算誤差には、丸め誤差と呼ばれるものがある。

たとえば、１０進数の0.1は２進数で書くと

という無限小数、循環小数になってしまう。

したがって二進数の小数点８桁で打ち切ると0.00011001（二進数）となる。この0.00011001（二進数）を十進数に直すと0.09765625（十進数）となるので、十進数で0.00234375という誤差が発生することになる。

このような計算機特有の誤差を丸め誤差という。

 

この他に桁落ちというものがある。

例えば、

という計算を考える。

いま仮に１０進数で８桁（浮動小数）しか精度を有さない計算機があるとする。

そして、これを使って上の計算を行おうとすると、

最初８桁あった有効精度が５桁に落ちてしまう。

このように、近い数同士の引き算をする場合、計算の精度が落ちてしまうことがある。

このような現象を桁落ちという。

しかし、次のように計算すると、

と精度が保たれる。

今のコンピュータ、とりわけ電卓は賢いので、どのように計算しても同じ計算結果を出すが・・・。

 

わかりやすいように１０進数を使って説明しているけれど、正確な議論をする場合、コンピュータは２進数で計算をしているので、２進数を使って議論をしなければならない。

この点は留意して欲しい。

 

この他に、たとえば、十進数で３桁しか精度を持たないコンピュータの場合、1000+1を計算すると、

となってしまう。

このように、大きな数（？）と小さな数（？）の足し算でも情報が落ちる情報落ちという現象が起きる。

整数３桁しか扱えないコンピュータなどあるものかと思うかもしれないけれど、かつて存在した８ビットのパソコンで扱える整数は、特別な処理を施さないかぎり、±の符合を有する整数の場合−128〜+127、符合を有さない整数の場合0〜255だ。

という計算をするためには、符号なしの整数の場合、最低でも10ビット必要でそれほど荒唐無稽な話ではないのである。

上の有効精度３桁というのは極端な例だが、Fortranで1000000+0.001を単精度で計算させると、次のような計算結果が得られる。

 

 

 

0.001という情報が落ちてしまっている。

つまり、情報落ちが発生している。

また、先に話したように、0.001は、一旦、２進数に変換されるので、これを１０進数に戻すときに1.00000005×10⁻³となっており、丸め誤差が発生していることもわかると思う。

 

このように、コンピュータで計算する場合、打ち切り誤差以外の様々な誤差が計算に混入する。そして、時に、計算途中で誤差が誤差を次々と生み出し、それが雪だるま式に膨れ上がり、この誤差の蓄積のために、とんでもない計算結果が得られるのであった。

 

タグ：数値解析

［境界要素法プログラムを設計する（基本のGauss掃き出し法）］

［境界要素法プログラムを設計する（基本のGauss掃き出し法）］

 

１．Gauss掃き出し法の標準構成

　前回の[外部手続き]の最後に、

    REM Gauss掃き出し法 **********************************************
    External Sub Gauss_0(A, z, n, Zero_Order, Return_Code)
      REM まだ何にもしない
    End Sub

という、よくわからんものがくっついてたはずです(^^)。今回はこの中身をつくろうという話です。Gauss掃き出し法は連立一次方程式を解く、標準解法です。なので別に境界要素法でなくても良いのですが(^^;)。

　Gauss_0の引数Aは、解くべき連立一次方程式の係数マトリックスを表す配列，zはGauss_0が呼ばれた時は連立一次方程式の既知ベクトル（配列）で、Gauss_0の処理が終わると解ベクトルになっています。nは未知数の数，Zero_Orderは後述するピボット選択で使用する0判定オーダーを表す整数，Return_Codeは計算終了状態をGauss_0の外部へ伝えるために用意した引数です。

 

　標準解法なので、数値計算業界（？）の業界標準みたいな構成があります。こういう構成です。

 

    External Sub Gauss_0(A, z, n, Zero_Order, Return_Code)
      External Sub Scaling
      External Sub Foward_Eliminate
      External Sub Back_Sub

      REM ｽｹｰﾘﾝｸﾞ
      Call Scaling(A, z, n, Return_Code)

      If Return_Code <> 0 Then
        Return
      End If

      REM 前進消去
      Call Foward_Eliminate(A, z, n, Zero_Order, Return_Code)

      If Return_Code <> 0 Then
        Return
      End If

      REM 後退代入
      Call Back_Sub(A, z, n)
    End Sub

 

　とりあえず、If Return_Code <> 0 ThenのIfブロックはおいといて、REM文の後のサブルーティン呼び出しに注目です。呼ばれるサブルーティンは、スケーリング(Scaling)，前進消去（Foward Elimination），後退代入（Back Substitution）を行います。このうち前進消去と後退代入の意味がわからない人は、

ネコ先生の記事を読もう！。

 

２．スケーリング

　スケーリングは桁落ち誤差対策です。連立一次方程式の係数マトリックスを表す配列Aは、ふつう倍精度浮動小数点実数で宣言されますが、倍精度浮動小数点実数型の公称有効数字は16桁と決まっています（国際規格）。以下は極端な例ですが、次の2元連立一次方程式を考えます。

　　

　上記を解くと、(2)－(1)で、

　　

　(2')より、

　　

なので、(3)を(1)に代入し、

　　

となります。(3)，(4)より厳密解は、ほとんどx＝2，y＝1である事がわかります。数値計算には必ず数値誤差が憑きもの（？(^^;)）です。なので例え数値誤差が憑いたとしても、x＝2，y＝1程度の近似解は出るようにしたい訳です。

　同じ事をコンピューターにやらせると、既に(2')の段階でコンピューター内部では-(1+10²⁰) → -10²⁰，1-10²⁰ → -10²⁰という現象が起きてしまいます。有効数字が16桁しかないからです。

　すると(3)でy＝1となり、それを(1)に代入するとx＝0となって全くお話になりません。何故ならx＝0，y＝1を(2)に代入すると、誤差がありとしても想定外の、-1＝1という不当過ぎるな結果になるからです。

　一方、

　　

としてもOKなはずです。(1')は(1)の両辺を10²⁰で割っただけのもの。ここで(2)－(1')をやるために、(1')のxの係数を1にしようとして、(1')を10²⁰倍したら馬鹿です。(1)と(2)の状態に戻るので。(1')－10^-20×(2)を試してみましょう。下から上を引いた方が考えやすいので行交換を行い、

　　

として(1')－10^-20×(2)を行うと、

　　

　(1")より、

　　

 

と(3)と同じものが得られますが、これを(2)に代入すれば、x＝2が得られます。

　だまされたような気持ちですか？(^^)。でもこのトリックは正しいんです。理由はなんでしょう？。連立方程式の「構造」に注目します。(1)と(1')を比べると、条件(1')ではxの寄与は非常に小さいとわかります。無視してもいいくらいです。無視するとy＝1です。それを(2)に代入しx＝2。(1)と(1')は同じ条件のはずなのに・・・。

 

　以上を強引にまとめると、こうなります。

　大きすぎる数値で桁落ちが起こるより、非常に小さい数値で桁落ちが起こった方がましだ！。

　そして桁落ちは必ず起こる！。

 

　理由は連立方程式の中で微小な係数しか持たない部分は、最初から無視しても構造的に被害小だからです。いいかえれば、ちゃんとした近似構造を解いた。

　しかし連立方程式の各条件は、何倍しても理屈の上では同じ条件。本当は寄与が少ない部分も、生の数値ではそう見えないかも知れない。よって寄与が少ないと見分けるためには、連立方程式の各行の数値桁数を揃えればいい、という話になります。そのために、連立方程式の各行を各行の絶対値最大成分の値で割り、各行の最大係数を1に揃える作業を、スケーリングと言います（←　一種のスケール規模問題です）。

　スケーリングを行うと、連立方程式の数値的安定性は格段に向上します。でも数値誤差は憑き物(^^;)。弊害もあります。ほとんど特異な連立方程式（解けちゃいけない）が、安定に解けてとんでもない答えが出る場合もあります。そして人間は馬鹿だから、その答えを信じます(^^;)。

 

　ところで(1')(2)を(2)(1')にしたような行交換を、機械は自動で出来るものでしょうか？。それは次のピボット選択と関連します。

　ところで桁落ち誤差の意味がわからない人は、ネコ先生の記事を読もう！。

 

 

 

３．ピボット選択

　ピボット選択は基本的に丸め誤差対策です。掃き出し中の連立一次方程式の係数配列は、概ね下図のようになります。灰色の部分は掃き出し完了ブロックで、成分は全て0です。Active行，Active列に囲われた部分が、これからs回目の掃き出しを行おうという、未処理ブロックになります。Active行，Active列の交点を、特にピボット（Pivot：枢軸）と言います。

 

　掃き出し処理は、Active列のs+1行～n行成分を0にすれば完了です。そのためにまず、Active行をピボットa_ssで割り、被処理行の頭の値a_{s+2 s}などをActive行にかけて、非処理行から引きます。「/」を使う時には、常に割る数が0でないかどうかのチェックが必要です。

　a_ssがたまたま0という事はあり得ます。でもActive列の中には非零成分がたぶんあります。そこでActive列のs行～n行を検索し、最初に見つかった0でないa_ksを持つ行をActive行（s行）と交換し、k行目をActive行とすればOKです。もしActive列に非零成分がなければ、特異行列です。答えがでちゃいけません。そこで処理を中止します。

 

 

　でも数値誤差は憑き物。本当は0なんだけど0に近いが0になってない数値はきっと、Active列にもごろごろしてます(^^;)。どれくらい小さければ0とみなすか？。要するに閾値Epsはどれくらいか？。これは問題のスケール規模で決まります。

　係数配列の数値は平均的に非常に大きい場合もあるし、逆もありえます。両方で同じEpsを使っては計算精度なんか保証されません。でもスケーリングしたのでした。各行の最大成分は初期状態で必ず1です。スケール規模は「1」です。「1」を基準にEpsを決定できます(^^)。ここで引数Zero_Order登場です。Active列の非零成分検索で、絶対値が　Eps＝1.0 * 10^(-Zero_Order)　未満のものは0とみなします。自分は経験的にZero_Order＝8を使います。

 

　でも非零という条件だけでいいのでしょうか？。係数配列は初期状態においてすら、すでに各数値には丸め誤差δが含まれます。コンピューターの数値には、有限の有効桁数しかないからです。コンピューターの中の方程式とは、要するに最初から近似方程式なんですよ。数値誤差は数値計算の憑き物！(^^;)。

　真の値がa_ij＋δならば、それをピボットa_ssで割る事により、丸め誤差の影響はδ/a_ssになります。小さすぎるa_ssで割ると、丸め誤差の影響は拡大します。もう一つ困るのはa_ssが微小だと、スケーリングのところでやった(1')を、わざわざ(1)に戻すような事まで起きる可能性がある事です。せっかくスケーリングしたのに。

 

　なので、Active列の絶対値最大成分を持つ行をActive行と交換します。これがピボット選択（行ピボット選択）です。あれっ、(1')と(2)の行交換が自動で出来ちゃった・・・(^^)。

 

　スケーリングとピボット選択をセットで使うと、概ね6桁くらいの有効数字を数値誤差の影響から救えます。10^6＝1,000,000なので、馬鹿にならない数字です。

　ところで丸め誤差の意味がわからない人は、ネコ先生の記事を読もう！。

 

REM ｽｹｰﾘﾝｸﾞ
External Sub Scaling(A, z, n, Return_Code)

  For i = 1 to n
    Ma = 0
    jj = 0

    REM 行の絶対値最大成分を検索
    For j = 1 to n
      h = Abs(A(i, j))

      If Ma < h Then
        Ma = h
        jj = j
      End If

    Next j

    REM 行が0だったら処理中止
    If Ma = 0 Then
      Return_Code = 1
      Return
    End if

    REM 絶対値最大成分で行を割り規格化
    h = A(i, jj)
    z(i) = z(i) / h

    For j = 1 to n
      A(i, j) = A(i, j) / h
    Next j

  Next i

  Return_Code = 0
End Sub

 

REM 前進消去
External Sub Foward_Eliminate(A, z, n, Zero_Order, Return_Code)

  REM 0判定値セット
  Eps = 10^(- Zero_Order)

  For i = 1 to n

    REM Pivot検索
    Ma = 0
    ii = i

    For k = i to n
      h = Abs(A(k, i))

      If Ma < h Then
        Ma = h
        ii = k
      End If

    Next k

    REM Pivot候補がEps未満だったら処理中止
    If Ma < Eps Then
      Return_Code = 2
      Return
    End if

    REM 行選択Pivot
    If ii <> i Then
      P1 = z(i)
      P2 = z(ii)
      z(i) = P2
      z(ii) = P1

      For j = i to n
        P1 = A(i, j)
        P2 = A(ii, j)
        A(i, j) = P2
        A(ii, j) = P1
      Next j

    REM 掃き出し処理
    h = A(i, i)  REM Pivotｾｯﾄ
    z(i) = z(i) / h  REM 既知ﾍﾞｸﾄﾙのActive行を規格化

    For j = i + 1 to n  REM 係数配列のActive行を規格化
      A(i, j) = A(i, j) / h
    Next j

    For k = i + 1 to n  REM 掃き出し
      h = A(k, i)  REM 非処理行の頭をｾｯﾄ

      If h <> 0 Then  REM Fillin Skip
        z(k) = z(k) - h * z(i)  REM 既知ﾍﾞｸﾄﾙの掃き出し

        For j = i + 1 to n  REM 係数配列の被処理行を掃き出し
          A(k, j) = A(k, j) - h * A(i, j)
        Next j

      End If

    Next k

  Next i

  Return_Code = 0
End Sub


REM 後退代入
External Sub Back_Sub(A, z, n)

  For i = n to 1 Skip -1
    h = 0

    For j = n to i + 1 Skip -1
      h = h - A(i, j) * z(j)
    Next j

    z(i) = h
  Next i

End Sub

 

(執筆　ddt²さん）

タグ：境界要素法 数値解析

極方程式で与えられた曲線に囲まれた部分の面積

極方程式で与えられた曲線に囲まれた部分の面積

 

曲線が極座標r、θを用いて

　　

で表されているとき、この方程式を曲線Cの極方程式という。

 

例えば、半径a（a>0）で点(a,0)を中心とする円

　　

は、

　　

と、極方程式で表すことができる。

 

また、デカルト直交座標における動点Pの座標を(x,y)とすると、

　　

であり、

　　

という関係がある。

 

 

定理

曲線Cがで表されていて、f(θ)が連続であるとする。このとき、曲線Cと半直線θ=α、θ=βで囲まれた部分の面積は

　　

である。

【略証】

　　

（証明終）

 

 

問題１　次の曲線（カージオイド）に囲まれている部分の面積を求めよ。

　　

【解】

この曲線はx軸に関して対称なので、上半分の面積S₁を求め、それを２倍すればよい（右図参照）。

したがって、

　　

したがって、この曲線によって囲まれる面積Sは

　　

（解答終）

 

問題２　次の曲線（レムニスケート）に囲まれている部分の面積を求めよ。

　　　　

【解】

この曲線はx軸、y軸に関して対称。だから、第１象限の面積を４倍すればよい。

　　

とおくと、

　　

は、

　　

rが恒等的に0のときは曲線ではなく、原点になるので、

　　

第１象限だけを考えればよいので、このとき、0≦θ≦π/2であり、また、（２）を満たさなければならないので、

　　

よって、

　　

したがって、

　　

（解答終）

 

 

問題３　次の曲線（３葉線）に囲まれる部分の面積を求めよ。

　　

【解】

とおくと、

求める面積は、第１象限でこの曲線によって囲まれる部分の面積S₁を３倍したもの。

ところで、0≦θ≦π/2で

　　

になるのは、

　　

したがって、

　　

よって、

　　

（解答終）

 

タグ：重積分 定積分 微分積分

テーラー級数、マクローリン級数を使って

テーラー級数、マクローリン級数を使って

 

まずは、テーラー級数を用いる次の問題を。

 

問題１

　　

をx=1のまわりのテーラー級数に展開せよ。

これを利用して、log1.1を小数３桁まで正確に求めよ。

【解】

　　

したがって、テーラー級数は

　　

または、剰余項を用いて

　　

log1.1を小数点以下３桁まで求めるには剰余項を

　　

にすればよいので

　　

したがって、n=3にとればよい。

よって、

　　

で、誤差は

　　

（解答終）

 

電卓で自然対数log1.1を求めると

　　

したがって、小数点以下３桁まで正確に求められており、誤差の評価も概ね正しいことがわかる。

 

 

問題２　マクローリン展開の最初の４項を用いて次の微分方程式

　　

の近似値を求めよ。

【解】

　　

の両辺を順次xで微分すると、

　　

したがって、

　　

よって、

　　

（解答終）

 

 

問題３　次の微分方程式を解け。

　　

【解】

この微分方程式はn=2のベルヌーイ形の微分方程式。

だから、

　　

とおくと、

　　

よって、

（解答終）

 

グラフを見ると、問題２で求めた近似解は、−0.5≦x≦0.5の範囲で厳密解とよく一致していることがわかる。

なお、

　　

タグ：数値解析 微分積分

差分法の基礎

差分法の基礎

 

f(x)を何回でも微分可能な関数とすると、テーラーの定理より、f(x)はx=aの近傍で

と展開することが可能である。

ここで、記号Oはランダウの（ビッグ）Oで

　　

ある。

 

したがって、n=1,2,3とすると、

　　

が成り立つ。

また、hを−hに置き換え、

　　

を得ることができる。

 

（３）式より、h≠0のとき、

　　

となり、このことからf'(a)を

　　

と近似したときの誤差はhのオーダー、であると予想できる。

 

このことは、拡張された平均値の定理

　　

より、h≠0のとき、

　　

となり、f''(x)は

　　

で連続だから

　　

となるMが存在し、

 　　

となり、

　　

と近似した誤差が｜h｜のオーダー程度であることが確かめられる。

f(x)がxの１次関数の時、

　　

なので、

　　

が成立すので、（６）は１次の精度である。

 

さらにa=0とし、hを変化させ、

と近似したときの誤差を求めてみると、右の図のようになる。

横軸にはh、縦軸に誤差をとり、対数グラフで結果を表している。

この直選の傾きが1であることから、（６）の近似式の誤差がhの１次オーダーであることが確かめられる。

このことを、

　　

と表すことにする。また、f'(a)のこの近似式を前進差分と呼ぶことにする。

 

また、f'(a)を求めるために、（３）と（３’）の辺々を引くと、

　　

が得られる。

（８）式で与えられるf'(a)の近似式を中心差分と呼び、この近似式の誤差はhの２次のオーダーであることと予想できる。

さらにa=0とし、hを変化させ、

　　

と近似したときの誤差を求めると右図のようになる。このときの直線（緑色の直線）の勾配は２であり、この近似式の誤差はhの２次オーダーであることがわかる。

f(x)がxの２次関数であるとき、

　　

となるので、

　　

が成立する。したがって、この近似式は２次の精度を持っている。

 

要するに、という記号は、n次の精度で、誤差の程度はのオーダーであり、hを1/10にすると、誤差はになるということを表していると考えることができる。

 

の近似式を求めるために、（５）と（５’）の辺々を加えると、

　　

したがって、

　　

という近似式は、２次の精度をもっており、誤差はh²程度ということになる。

 

この直線の傾きは２であり、誤差がhの２次のオーダーであることが確かめられる。

タグ：微分積分 微分 数値解析

連立方程式の解法３ 反復法１

連立方程式の解法３　反復法１

 

ヤコビ法

 

次の３元連立方程式があるとする。

　　

①をx、②をy、③をzについて解くと次のようになる。

　　

ここで、として上の式の右辺に代入すると、

　　

となる。

このを用いて

さらに、この値を用いて

　　

このように繰り返し計算すると、

と、やがて、３元連立方程式の解である(x,y,z)=(1,2,3)に収束してゆく。

 

このように反復計算をして連立方程式の解を求める方法をヤコビ法という。

 

より一般の

　　

つまり、

　　

の場合、次のように計算すればよい。

 

 

STEP1　１行からn行まで

　　

を計算する。

 

STEP2　i=1〜nまでのすべてのをの値に更新する

 

STEP3　所定の精度に達していなければSTEP1に戻る。達していたら、計算を打ち切る。

 

 

いま述べたヤコビ法は、直接、連立方程式の解を求める掃き出し法やガウス消去法とは異なり、連立方程式が唯一の解を持っていたとしても、計算が収束しなかったり、発散する場合がある。

優対角条件、つまり、

　　

などの条件が成立しないとき、ヤコビ法などの反復法は振動したり、発散する。

 

このことは

　　

の②と③を入れ換え

　　

をヤコビ法で解かせると、このことが確かめられる。

 

C言語によるサンプルプログラムを次に示す。

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


void Jacovi(double a[][N], double x[], double b[], int n) {
    double eps = 1.e-6, sum=0, res=0;
    double xn[N];
    int i, j, k;
    int limit = 100;
   
    printf("反復回数 x(1)      x(2)      x(3)      残差\n");
    for (k = 1 ; k <= limit; k++) {
        res = 0;
        for (i=0; i < n; i++) {
            sum = b[i];
            for (j=0; j < n; j++) {
                sum = sum - a[i][j]*x[j];
            }
            res = fabs(sum);
            xn[i]= x[i]+sum/a[i][i];
        }
        if (k%5==0) { // ５回に１度、計算途中の結果を表示
            printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, xn[0],xn[1], xn[2], res);
        }
        for (i=0; i<n; i++) {
            x[i]=xn[i];
        }
        if (eps>res) break;
    }
}

main() {
    int i, n;
    double a[N][N] = {{3,1,1}, {1,3,1}, {1,1,3}};
    double x[N]={0.};
    double b[N]= {8., 10., 12.};
   
    n=3;
   
    Jacovi(a,x,b,n);

    for (i = 0; i < n; i++) {
        printf("x[%d] = %f\n", i+1,x[i]);
    }
}

タグ：数値解析

連立方程式の解法３　反復法１

 

ヤコビ法

 

次の３元連立方程式があるとする。

　　

①をx、②をy、③をzについて解くと次のようになる。

　　

ここで、として上の式の右辺に代入すると、

　　

となる。

このを用いて

さらに、この値を用いて

　　

このように繰り返し計算すると、

と、やがて、３元連立方程式の解である(x,y,z)=(1,2,3)に収束してゆく。

 

このように反復計算をして連立方程式の解を求める方法をヤコビ法という。

 

より一般の

　　

つまり、

　　

の場合、次のように計算すればよい。

 

 

STEP1　１行からn行まで

　　

を計算する。

 

STEP2　i=1〜nまでのすべてのをの値に更新する

 

STEP3　所定の精度に達していなければSTEP1に戻る。達していたら、計算を打ち切る。

 

 

いま述べたヤコビ法は、直接、連立方程式の解を求める掃き出し法やガウス消去法とは異なり、連立方程式が唯一の解を持っていたとしても、計算が収束しなかったり、発散する場合がある。

優対角条件、つまり、

　　

などの条件が成立しないとき、ヤコビ法などの反復法は振動したり、発散する。

 

このことは

　　

の②と③を入れ換え

　　

をヤコビ法で解かせると、このことが確かめられる。

 

C言語によるサンプルプログラムを次に示す。

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


void Jacovi(double a[][N], double x[], double b[], int n) {
    double eps = 1.e-6, sum=0, res=0;
    double xn[N];
    int i, j, k;
    int limit = 100;
   
    printf("反復回数 x(1)      x(2)      x(3)      残差\n");
    for (k = 1 ; k <= limit; k++) {
        res = 0;
        for (i=0; i < n; i++) {
            sum = b[i];
            for (j=0; j < n; j++) {
                sum = sum - a[i][j]*x[j];
            }
            res = fabs(sum);
            xn[i]= x[i]+sum/a[i][i];
        }
        if (k%5==0) { // ５回に１度、計算途中の結果を表示
            printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, xn[0],xn[1], xn[2], res);
        }
        for (i=0; i<n; i++) {
            x[i]=xn[i];
        }
        if (eps>res) break;
    }
}

main() {
    int i, n;
    double a[N][N] = {{3,1,1}, {1,3,1}, {1,1,3}};
    double x[N]={0.};
    double b[N]= {8., 10., 12.};
   
    n=3;
   
    Jacovi(a,x,b,n);

    for (i = 0; i < n; i++) {
        printf("x[%d] = %f\n", i+1,x[i]);
    }
}

ネムネコ、表計算ソフトで、平行平板間の流れ（平板ポアズイユ流れ）を解く！！

平行平板間の流れ

 

平行平板間の流れに関する、次の微分方程式を、表計算ソフトを使って解いてみたにゃ。

　　

境界条件は

　　

これがその結果。

 

 

 

 

ΔY=0.1という非常に粗い分割でも、相対誤差約１％で解けてしまうだケロ。

すごいもんだろ。

 

ネットにスプレッドシートを公開したんで、興味のあるヒトはアクセスするといいにゃ。

https://docs.google.com/spreadsheets/d/e/2PACX-1vStliIRWmWqLXymnL_LCGJFXDqbN5LGMPESBKVfojhSfHN5vcU1IqNvyip0sIT8Zm5YOfSupub1ii8i/pubhtml

エクセル版は↓
https://docs.google.com/spreadsheets/d/e/2PACX-1vStliIRWmWqLXymnL_LCGJFXDqbN5LGMPESBKVfojhSfHN5vcU1IqNvyip0sIT8Zm5YOfSupub1ii8i/pub?output=xlsx

LibreOffice calc版は↓

https://docs.google.com/spreadsheets/d/e/2PACX-1vStliIRWmWqLXymnL_LCGJFXDqbN5LGMPESBKVfojhSfHN5vcU1IqNvyip0sIT8Zm5YOfSupub1ii8i/pub?output=ods

だにゃ。

平行平板間の流れ（平面ポアズイユ流れ）

平行平板間の流れ（平面ポアズイユ流れ）

 

右図のような２つの無限平行板間の流れは、次の方程式で表される。

　　

ここで、pは圧力、uはx軸に平行な流速、さらに、Qは（体積）流量、μは粘性係数である。

なお、（１）の境界条件は

　　

 

（１）は、流体力学の基礎方程式であるナビエ・ストークス方程式

　　

の最も簡単なものの一つで、ナビエ・ストークスの解析解（厳密解）が得られるものの一つで、平均流速

　　

を用いると、解は

　　

である。

 

u、yを

　　

と無次元化すると、（５）式は

　　

となり、UはY=1/2のときに

　　

という最大値をとる。

さらに、pとxを次のように無次元化すると

　　

（１）、（２）式は

　　

ここで、Reはレイノルズ数と言われる、流体力学で流れを支配する、重要な無次元数で

　　

である。

 

実は、このReが大きくなると、ナビエストークス方程式の非線形項の非線形性が強くなり、流れは時間、空間的に不規則――不規則って何だ(・・?――な乱流という、数学・物理学的に厄介なものになるのだけれど・・・。

 

それで、今回、やりたいのは、差分法を用いて、（１０）と（１１）を連立方程式に直し、これを数値的に解こうということ。

 

（７）式を見るとわかるけれど、この流れの速度分布は、Y=1/2を軸とする放物線で表されており、Y=1/2で対称である。

したがって、

　　

だから、

（１０）式の境界条件は、

　　

ではなく、

　　

としてもよい。

 

なぜ、境界条件を、（１４）ではなく、（１５）とするかというと、計算領域を0≦Y≦1から0≦Y≦1/2と小さくし、計算量をできるだけ少なくしたいから。

 

[0,1/2]をn等分し、

　　

とおき、

　　

とし、この点における無次元流速Uをとする。

差分法を用いると

　　

となるので、（１０）式は

　　

となり、

　　

ここで、

　　

とおけば、（１９）式は

　　

になる。

Y=0の境界条件U=0より

　　

問題は、Y=1/2における境界条件

　　

をどのように扱うかだけれど、中心差分を用いると、

　　

となるので、i=nのとき

　　

とすればよいだろう。

 

（１１）式も対称性を利用すると、

　　

となり、左辺の積分を台形公式で近似すれば、

　　

となる。

よって、

　　

というn+1元の連立方程式が得られる。

未知数は、のn個とdP/dXの１個で、合せてn+1個だから、連立方程式（２５）を解くことによって、この未知数を全て求めることができる。

（２５）は行列で書き換えると、

　　

したがって、係数行列を

　　

にすることにより、ガウス消去法や掃き出し法を用いて連立方程式（２６）を解くことができる！！

 

（２６）の左辺の行列の成分の圧倒的多数は０の疎行列なので、単純なガウス消去法などの直接法は無駄な計算が多く、効率的ではないのだけれど・・・。

 

n=10の場合の計算結果は次のとおり。

 

 

速度U、圧力勾配dP/dXの計算誤差は約0.25％に収まっていることがわかる。

 

実は、ここまでのプログラムを作ると、それをもとに、より複雑なNS方程式の特殊形

　　

で表される無限平板間の流れを数知的に解くこともできる！！

 

なのですが、上の方程式は非線形な上に、計算点が少なくとも数百〜数千程度になるので、ガウス消去法などの直接法で、上の連立微分方程式から得られる連立方程式を解くなど狂気の沙汰！！

 

連立方程式

　　

を効率よく解く数値計算法、つまり、アルゴリズムを作らないといけない。

 

ガウス消去法の計算量はn³程度になるという話を前にしたけれど、ガウス消去法をもとに、計算量をnオーダーに留める計算手法、アルゴリズムを作ることは可能！！

我はと思うヒトは、この計算手法を考えてみるといいにゃ。

それさえできれば、

　　

を数値的に解くプログラムを作ることだって可能になる！！

タグ：数値解析

最小２乗法による近似式の決定

最小２乗法による近似式の決定

 

§１　最小２乗法による近似多項式の求め方

 

n個のデータがあるとする。

このデータがm次の多項式（m≦n+1）

　　

に沿っているとし、

　　

が最小になるようにm次多項式（１）の係数を定める方法を最小２乗法という。

最小２乗法では、

　　

とおいたとき、（２）式が最小になるようにを定めるので、

　　

にならなければならない。

したがって、

　　

でなければならない。

ところで、

　　

となるので、（４）式は

　　

ここで、

　　

とおくと、（５）式は

　　

となる。

したがって、この連立方程式（７）を解くことによって、多項式（１）の係数を定めることができる。

 

特に、m=1のとき、すなわち、

　　

のとき、（７）は

　　

ここで、

　　

となるので、（８）は

　　

これを解くと、

　　

ここで、

　　

したがって、統計学の回帰直線と一致する。

 

§２　プログラム

 

n個のデータをもとに、最小２乗法を用いてm次の近似多項式を求める連立方程式は（６）、（７）式から与えられる。

（６）、（７）から与えられる連立方程式をガウス消去法を用いて解くことにすると、例えば、次のようなプログラムを作ることができる。

 

real x(20),y(20) !xとyのデータは２０個まで
real c(0:10) ! 近似多項式の係数（１０次まで）
real y_c(20)
data x/20*0/ ! x(i)の初期化
data y/20*0/ ! y(i)の初期化

write(6,*) "データの個数nを入力"
read(5,*) n

write(6,*) "データの入力"
read(5,*) (x(i), i=1,n) ! xのデータをn個入力
read(5,*) (y(i), i=1,n) ! yのデータをn個入力

write(6,*) "近似多項式の次数mを入力"
read(5,*) m

call lsm(x,y,n,m,c)

write(6,*)
! 近似式の係数の表示
write(6,*) "近似多項式の係数"
do i=0, m
    write(6,100) 'c(',i, ')=', c(i)
end do

write(6,*)
! 計算結果
write(6,*) "  x(i)       y(i)        計算値"
do i= 1, n
    s = 0
    do j=0, m
        s = s + c(j)*x(i)**j
    end do
    y_c(i)=s
    write(6,110) x(i), y(i), s
end do

open(1,file='lsm.dat', status='replace')
do i=1,n
    write(1,120) x(i),y(i),y_c(i)
end do

100 format(a,i2,a3,f10.7)
110 format(f10.7,2x,f10.7,2x,f10.7)
120 format(f10.7,x,f10.7,x,f10.7)
end

! ここから下は変えないほうがよい(^^)
! 最小二乗法
subroutine lsm(x,y,n,m,c)
real a(0:50,0:51)
real x(20), y(20), c(0:10)
real s(0:20), t(0:10)

! 係数の初期化
do i=0, 2*m
    s(i) = 0
end do
do i=0, m
    t(i)=0;
end do

! 係数を計算
do i=1, n
    do j=0, 2*m
        s(j) = s(j)+x(i)**j
    end do
    do j=0, m
        t(j)= t(j)+x(i)**j * y(i)
    end do
end do

! 係数行列に値をセット
do i= 0, m
    do j=0, m
        a(i,j) = s(i+j)
    end do
    a(i,m+1) = t(i)
end do

! ガウス消去法で連立方程式を解く
call Gauss(a,m)

do i=0, m
    c(i) = a(i,m+1)
end do

end

subroutine Gauss(a,n) ! ガウス消去法
real a(0:50,0:51)

! foward marching process
do k=0, n-1
! ピボット選択
    i_max = k
    do j=k+1, n
        if (abs(a(k,i_max)).lt.abs(a(k,j))) i_max = j
    end do
    if (i_max.ne.k) then
        do j=k,n+1
            w=a(i_max, j)
            a(i_max, j) = a(k,j)
            a(k,j) = w
        end do
        det=-det
    end if
    if (abs(a(k,k)).lt.eps) then
        det =0.
        return
    end if
! ピボット選択終わり
    do i=k+1, n
        t=a(i,k)/a(k,k)
        do j=k+1, n+1
            a(i,j)=a(i,j)-t*a(k,j)
        end do
    end do
end do

! backward marching process
do i=n,0, -1
    d = a(i,n+1)
    do j= n, i+1, -1
        d=d- a(i,j)*a(j,n+1)
    end do
    a(i,n+1)=d/a(i,i)
end do
end

 

 たとえば、次の３個のデータ

(1,1) , (2,5), (3,8)

をもとに２次の近似多項式を求めると、次のようになる。

 

 

さらに、次の７個のデータ

(−3,5),(−2,−2),(−1,−3),(0,−1),(1,1)、(2,4),(3,5)

をもとに５次の近似多項式を求めると、次のようになる。

 

　

タグ：数値解析