So-net無料ブログ作成
ネコ騙し数学 ブログトップ
前の10件 | 次の10件

[境界要素法プログラムを設計する(境界積分サブルーティン:実装編)] [ネコ騙し数学]

[境界要素法プログラムを設計する(境界積分サブルーティン:実装編)]

1.積分パラメータ関数の実装

 

bem14-fig-001.png

 

 前々回、メインの中に与えた境界積分サブルーティン(関数)は、以下でした。

 

    REM  (x0,y0),(x1,y1),(x2,y2)から、要素kの境界積分値を計算
    REM  ψj1に関する結果をB1,qj1に関する結果をH1で表す
    REM  ψj2に関する結果をB2,qj2に関する結果をH2で表す

    Lk = Length(x1,y1,x2,y2)
    Beta_k = Angle(x2 - x1, y2 - y1, Lk)
    r1 = Length(x0,y0,x1,y1)
    r2 = Length(x0,y0,x2,y2)
    Gamm_1 = Angle(x1 - x0, y1 - y0, r1)
    Gamm_2 = Angle(x2 - x0, y2 - y0, r2)
    h = Normal_Length(x2 - x1, y2 - y1, x0,y0)
    s = Normal_Foot(x2 - x1, y2 - y1, x0,y0, h)

  B1 = B_Inte_B1(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
  B2 = B_Inte_B2(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
  H1 = B_Inte_H1(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)
  H2 = B_Inte_H2(Lk, r1, r2, Beta_k, Gamm_1, Gamm_2, h, s)

 B1B2H1H2の計算前に必要な積分パラメータを全て揃えます。そこには、B_Inte_?の中身として[内点方程式の積分公式(ラプラス型)]と[内点方程式から境界方程式へ(ラプラス型)]で導出した式を単純に並べたいという意図があります。なのでまず、積分パラメータを与える関数を検討します。

 

 最初にLength関数ですが、明らかに次でOKです。

 

    External Function Length(x1,y1,x2,y2)
      Length = Sqr((x2 - x1)^2 + (y2 - y1)^2)
    End Function

 前々回に予告したように、関数とサブルーティンはローカル変数を持てる、外部手続きで統一します。Normal_LengthNormal_Footの計算に疑問の余地はないのですが、ちょっと長そうなので後にします。

 次にBeta_kの扱いですが、前回の検討からGamm_1Gamm_2の特殊ケースとして扱って良いとわかりました。という訳でGamm_1Gamm_2の計算です。これらは単独に計算しては駄目で、Gamm_2Gamm_1の形で計算する必要がありました。それをGam_12で表します。

 

 まず特異点(x0y0)が要素節点(x1y1)(x2y2)に一致しないケースから考え始めると、わかりやすいと思います。以下、前回を参照しつつ・・・。

 

 計算開始の前提として、図-1r1r2および(x0y0)(x1y1)(x2y2)は引数渡しされるとします。また特異点と要素節点の一致/不一致判定のために、0距離判定値Epsも必要です。

 

    External Function Angle(x0, y0, x1,y1,x2,y2, r1, r2, Eps)

      If r1 < Eps Then
        x01 = 1
        y01 = 0
        x02 = x2 - x1
        y02 = y2 - y1
        r1 = 1
      ElseIf r2 < Eps Then
        x01 = x1 - x2
        y01 = y1 - y2
        x02 = 1
        y02 = 0
        r2 = 1
      Else
        x01 = x1 - x0
        y01 = y1 - y0
        x02 = x2 - x0
        y02 = y2 - y0
      End If

      Gum_12 = Acos((x01 * x02 + y01 * y02) / r1 / r2)
      Sign = Sgn(x01 * y02 - x02 * y01)
      Gum_12 = Gum_12 * Sign

      Angle = Gum_12
    End Function

 ElseEnd Ifの設定が基本です。

 r1 < Epsr10)の時は、(x01y01)(10)とするのでした。この時(x0y0)(x1y1)に一致するとみなせるので、(x02y02)(x2 - x1y2 - y1)です。そしてr11とします。

 r2 < Epsr20)の時は、(x02y02)(10)とするのでした。この時(x0y0)(x2y2)に一致するとみなせるので、(x01y01)(x1 - x2y1 - y2)です。そしてr21とします。

 

 以上より、積分パラメータを与える関数呼び出しは、

 

    Lk = Length(x1,y1,x2,y2)
    r1 = Length(x0,y0,x1,y1)
    r2 = Length(x0,y0,x2,y2)
    Gam_12 = Angle(x0, y0, x1,y1,x2,y2, r1, r2, Eps)
    h = Normal_Length(x0, y0, x1,y1,x2,y2)
    s = Normal_Foot(x0, y0, x1,y1,x2,y2)

になります(Normal_LengthNormal_Footの引数が前回のままでは不味いと気付いた(^^;))。

 

 Normal_Lengthの実装は以下です。Normal_Lengthは図-1の垂線距離hを求めます。

a = x2 - x1

b = y2 - y1

とすると、(x1y1)(x2y2)を通る直線は、yb/a*xに平行です。a0の場合も考慮して、これをbxay0の形にし、(x1y1)を通る事を考慮すると要素kの直線の方程式は、b(xx1)a(yy1)0になります(なんか受験を思い出すなぁ~(^^))。直線の式からその垂線ベクトルの方位は、明らかに(b,-a)です。従って点(x0y0)からの垂線は、tを任意の実数として、

  

と書けます。(1)b(xx1)a(yy1)0に載る条件は、(1)を代入し、

  

 (2)tについて解くと、

  

 従ってhは、t(3)の値の時、

  

 

 tを利用すれば垂線の足もわかるので、図-1sも出せます。垂線の足の座標は、t(3)の値にした(1)です。よってsは、

  

から、

  

です。

 こうなるとNormal_LengthNormal_Footで同じパラメータtを利用したいですよね。こうしましょう。

 


    t = Normal_Param(x0, y0, x1,y1,x2,y2)
    h = Normal_Length(x1,y1,x2,y2, t)
    s = Normal_Foot(x0, y0, x1,y1,x2,y2, t)

    External Function Normal_Param(x0, y0, x1,y1,x2,y2)
      a = x2 - x1
      b = y2 - y1
      t = (b * (x1 - x0) - a * (y1 - y0)) / (a^2 + b^2)

      Normal_Param = t
    End Function

    External Function Normal_Length(x1,y1,x2,y2, t)
      a = x2 - x1
      b = y2 - y1
      h = Abs(t) * Sqr(a^2 + b^2)

      Normal_Length = h
    End Function

    External Function Normal_Foot(x0, y0, x1,y1,x2,y2, t)
      a = x2 - x1
      b = y2 - y1
      s = Sqr((t * b + x0 - x1)^2 + (-t * a + y0 - y1)^2)

      Normal_Foot = s
    End Function


2.境界積分関数の実装

 以上で境界積分に必要な積分パラメータは全て揃いました。

 

  B1 = B_Inte_B1(Lk, r1, r2, Gam_12, h, s, Eps)
  B2 = B_Inte_B2(Lk, r1, r2, Gam_12, h, s, Eps)
  H1 = B_Inte_H1(Lk, r1, r2, Gam_12, h, s, Eps)
  H2 = B_Inte_H2(Lk, r1, r2, Gam_12, h, s, Eps)


で良いはずです。引数の最後にEpsが入るのは、これらでも特異点と要素節点の一致/不一致判定が必要だからです。

 

 後はB_Inte_?の中に、[内点方程式の積分公式(ラプラス型)]と[内点方程式から境界方程式へ(ラプラス型)]で導出した式を並べるだけですが、関数の内部構造はExternal Function Angleと同一構造になるとわかります。なのでIfブロックの形を決め、「関数名 = 値」の行を書いてから、計算式を並べる事をお奨めします。



    External Function B_Inte_B1(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        B1 = Gam_12 / 2 /Pai
      ElseIf r2 < Eps Then
        B1 = 0
      Else
        B1 = Gum_12 - 1 / Lk * (s * Gum_12 + h * (Log(r2) - Log(r1)))
        B1 = B1 / 2 /Pai
      End If

      B_Inte_B1 = B1
    End Function

    External Function B_Inte_B2(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        B2 = 0
      ElseIf r2 < Eps Then
        B2 = Gam_12 / 2 /Pai
      Else
        B2 = 1 / Lk * (s * Gum_12 + h * (Log(r2) - Log(r1)))
        B2 = B1 / 2 /Pai
      End If

      B_Inte_B2 = B2
    End Function

    External Function B_Inte_H1(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        H1 = Lk / 4 / Pai * (Log(Lk) - 3 / 2)
      ElseIf r2 < Eps Then
        H1 = Lk / 4 / Pai * (Log(Lk) - 1 / 2)
      Else
        H1 = -1 / 4 / Pai / Lk * (r2^2 * Log(r2) - r1^2 * Log(r1) - 1 / 2 *((Lk - s)^2 - s^2))
        H1 = H1 + 1 / 2 / Pai * (1 - s / Lk) * (h * Gum_12 + (Lk - s) * Log(r2) + s * Log(r1) - Lk)
      End If

      B_Inte_H1 = H1
    End Function

    External Function B_Inte_H2(Lk, r1, r2, Gam_12, h, s, Eps)

      If r1 < Eps Then
        H2 = Lk / 4 / Pai * (Log(Lk) - 1 / 2)
      ElseIf r2 < Eps Then
        H2 = Lk / 4 / Pai * (Log(Lk) - 3 / 2)
      Else
        H2 = 1 / 4 / Pai / Lk * (r2^2 * Log(r2) - r1^2 * Log(r1) - 1 / 2 *((Lk - s)^2 - s^2))
        H2 = H2 + 1 / 2 / Pai * s / Lk * (h * Gum_12 + (Lk - s) * Log(r2) + s * Log(r1) - Lk)
      End If

      B_Inte_H2 = H2
    End Function

 以上の結果は、[計算部][外部手続き]のところへは、まだはめ込みません。[計算部]では、まだやる事があるからです(もぉ~、境界要素法って面倒くさいなぁ~(^^;))。


nice!(0)  コメント(0) 

2階常微分方程式の境界値問題を選点法で解く2 [ネコ騙し数学]

2階常微分方程式の境界値問題を選点法で解く2

 

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

【答】

 

前回、この微分方程式の解を

と近似し、残差

をとし、重み関数wにディラックのデルタ関数

をとり、

となるようにを定める有限要素法の選点法を用い、選点を1/2とすることにより、

という近似解を得た。

 

しかし、こうして得た近似解の誤差が大きかったので、より精度のよい近似解を求める方法を考える。

 

そこで、(1)式を

とおくと、

となる。

このとき、

となるので、残差は

となる。

未知数がa₁a₂の2つなので、これを定めるために選点をx=1/4x=1/2に選び、その位置での残差を0とすると、

fem2-001.png

(10)と(11)から

fem2-002.png

という連立方程式が得られ、これを解くと、

となり、

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

 

fem2-graph-001.png計算結果を右図に示す。

(12)は(6)よりも厳密解(1)からの乖離が小さくなり、(1)の挙動をよく捉えていることが分かる。

 

(9)式の

とおくと、

と表すことができる。このψ₁ψ₂試行関数(Trial Functionと呼ぶ。

さらに、

とおくと、微分方程式は

と書くことができ、残差は

で表すことができる。

より一般に書くと、

これに重み関数をかけ、

となるように、この連立方程式(15)から、未知数を定め、微分方程式(13)の近似解を有限要素法重み付き残差法という。

この重み関数のディラックのデルタ関数

を用いるものを選点法といい、デルタ関数の性質から

になる。

n=2のとき、未知数はa₁a₂の2つだから、点をx₁x₂の2点を選び、連立方程式

を解くことによって、a₁a₂を定めることができる。

 

これが計算の仕組みというわけ。

 


タグ:数値解析
nice!(0)  コメント(0) 

[境界要素法プログラムを設計する(境界積分サブルーティン:具体調査編)] [ネコ騙し数学]

[境界要素法プログラムを設計する(境界積分サブルーティン:具体調査編)]

 

1.最後はボトムアップになる

 

bem14-fig-001.png

 

 例は1個で十分です。

  

 

 式(1)は、[内点方程式の積分公式(ラプラス型)]であげた、境界要素kの節点jψjの係数を与える積分結果です(ずいぶん前のような気がしてますが(^^;))。特異点が要素節点に一致しない場合、境界方程式でもこの形の計算は必ずします。というかこういう場合の方が、圧倒的に多いです。

 式中のLksγ2γ1hr2r1が図-1から定められる積分パラメータです。このうちLkshr2r1は単純に計算できます。例えばr1なら明らかに、

 

  

 

ですし、shはちょっと手がかかりますが、基本は(x0y0)(x1y1)(x2y2)で係数が決まる22列の連立方程式を解くだけです。問題はγ2γ1の与え方です。これは一見簡単そうに見えます。例えばγ1なら、

 

  

でいいんでないの?と。まず汎用言語の数学関数のAtntan-1の事)は、値域がふつう-π/2π/2です。なんか鬱陶しそうだから、x1x0y1y0の符号に注意して場合分けすれば、容易に0をカバーできます。これでγ2γ10は容易に計算できそうです。でも図-2のような場合はどうしますか?。

 

bem14-fig-002.png

 

 Atn0をカバーできるように改造したとしても、図-2では明らかにγ2γ1γ2γ10です。しかし線積分(境界積分)は常に左回りなので、図中の角度矢印で示したように(左回り)、γ2γ10でなければならず、しかも0で表した|γ2γ1|は本当の(γ2γ1)と値が違います(大きい)。

 では境界要素がx軸方向を横切る時だけ、引いて符号を反転させればOKなのか?というと、一般的には図-2への対処も必要です。

 

bem14-fig-003.png

 

 この場合は角度矢印が示すように(右回り)、γ2γ10が正解です。図-3では0の角度定義でたまたまγ2γ1なので、そのままやって正解になりますが、図-23の複合ケースは?。・・・もう鬱陶しくてやってられません(^^;)

 

 つまりγ1γ2は(βkも)単独で計算しては駄目なんですよ。あくまでγ1γ2の角度差分を局所情報に基づいて計算する必要があります。こういう事は、普通の積分では起こりません。普通は図-1のように一周したら値は元に戻って問題ないはずです。これは積分結果がtan-1になるような正則でない関数を積分せざる得ない、特異積分を使っている境界要素法の特徴の反映です(複素積分を知ってる人なら、イメージできるはず)。そういう事情から、前回メインで大雑把に決めた関数の引数並びも、角度差分を与えるように変更する必要があります。

 これはボトムアップです。このように最後はボトムアップになる事は多いです。トップダウンだけでは問題の個別状況の全ては把握できないからです。でも最後までボトムアップは我慢した方が、たいてい上手く行きます。

 

 

2.角度差分を計算するアルゴリズムを考える

 一つだけ確かなのは、|γ2γ1|π以下だ、という点です。|γ2γ1|πに成り得るのは、特異点が境界要素の中点上などにある場合なので、πを越える事はありません。そうすると格好の数学関数があります。Acoscos-1)です。Acosの値域は、汎用言語においてたいてい0πです。内積の定義から、

 

  bem14-001.png

 

として、

  

なので

  

と求められます。|γ2γ1|の符号は要するに、図-23の角度矢印の向きがわかれば良いので、ベクトル(x01y01)(x02y02)の外積を利用する手があります。(x01y01)γ1方向,(x02y02)γ2方向のベクトルなので。(x01y01)(x02y02)3次元化して外積を取ると、

 

  bem14-002.png

 

になります。(7)右辺のz成分が+なら左回りで、γ2γ1|γ2γ1|(7)右辺のz成分が-なら右回りで、γ2γ1=-|γ2γ1|と判定できます。ただし式(6)(7)は小さすぎるr1r2に対しては、結果が不安定になるのは目に見えています。本当にr1r20なら、division by zeroエラーになります。なりますが、特異点が要素節点に一致する場合も、もともと扱うのでした。だとすればそういう特殊ケースは、特異点が要素節点に一致するケースに含めれば十分です。

 

3.特異点が要素節点に一致するケース

 

bem14-fig-004.png

 

 ここでもγ2γ1の扱いが問題になります。[内点方程式から境界方程式へ(ラプラス型)]でFree Termの項は、次のように計算されました。

 

  

 

 式(8)上段の1/2π×()内の1項目が要素k+1γ2γ1すなわち、図-4γ3γ(j+1)です。2項目は要素kγ2γ1すなわち、図-4γ(j+1) γ1です。ただしε→0の極限で考えてます。

 でも待って下さい。ε→0の場合(特異点が要素節点に一致する場合)、どうやってγ(j+1) を与えれば良いんでしょう?。だって数値計算上はε0なんですから。そしてすぐに気づくのは、特異点が要素節点に一致するケースでは、「個々の要素の境界積分値は、特異点の要素節点への近づき方に依存する!」という驚愕の事態が発覚します。もろな特異積分!。境界積分は結局、可積分でない?(^^;)

 ・・・とけっこう猛烈に焦っちゃう訳ですが、式(8)下段では不定なγ(j+1)が引き算でキャンセルされて、「境界積分全体では近づき方に依存せず、値は一意に定まる」事になります。境界要素法は、こういう危ういところを持っています。ブレビアさんがなりふり構わず、コーシーの主値積分を持ち込んだ気持ちもわかりますよ(^^;)

 

 つまり足せばγ(j+1)は消えるのですから、逆に言えばなんでも良い訳です。例えば要素k(x02y02)(10)とし、要素k+1でも(x01y01)(10)にするとか。共有節点でγの方法が揃ってればOKなはず。どうしてかと言うと特異積分値は、局所情報に基づいて計算する必要があるから。

 特異点が要素節点に一致するケースでは、その点でγ1γ2(10)方向に取ると決定します。後で後悔するかも知れませんが(^^;)

 

 あと考えるべき事は、r1またはr2がどれくらい0に近かったら、0とみなすかです。これは【桁落ち誤差と丸め誤差はプログラマーの敵!】のところで書いたように、

  http://nekodamashi-math.blog.so-net.ne.jp/archive/c2306081239-2

 

「問題のスケール規模」で決めるのが妥当と考えられます。r1r2は特異点と要素節点の距離です。これに対応するスケール規模は明らかに、解析領域の大きさです。0判定Epsはどこで与えれば良いでしょう?。

 手がかりはたいてい「実行者」と「妥当な実行タイミング」にあります。「実行者」の第一候補は[入力部]です。[入力部]Roll(役割)は、解析領域という外部情報をセットし[計算部]へ渡す事です。Epsは解析領域の大きさで決まるので、これは外部情報のセットです(そして渡す)。

 次に実行タイミングですが、Epsのセットは[計算部]の実行前に、事前に一回だけで十分です。[計算部]の実行前には、[入力部]の実行があります。「実行者」と「妥当な実行タイミング」が理想的に一致したので、[入力部]の中でEpsのセットを行います。セットの場所は、全ての入力が完了した後の方がなんか安心じゃないですか?。

 

 [メモ欄][宣言部]はできれば一対一に対応させたいので、[メモ欄]の最後に「REM Eps0判定,Eps00オーダー」を追加し、[宣言部]の最後には「REM Eps0判定」と「Eps0 = 10^(-6) REM 0オーダー」のコードを追加します。コードの心は、解析領域の大きさの1/100万より小さい距離は0とみなすです。そして[入力部]の最後に、

 

    REM 全ての入力完了
    REM 解析領域のx方向幅とy方向幅を取得
    x_Min = 10^10    REM x座標最小値初期化
    x_Max = -10^10   REM x座標最大値初期化
    y_Min = 10^10    REM y座標最小値初期化
    y_Max = -10^10   REM y座標最大値初期化

    For j = 1 to BN_Count
      h = BN(j, 1)  REM j節点のx座標

      If h < x_Min Then
        x_Min = h
      End If

      If x_Max < h Then
        x_Max = h
      End If


      h = BN(j, 2)  REM j節点のy座標

      If h < y_Min Then
        y_Min = h
      End If

      If y_Max < h Then
        y_Max = h
      End If

    Next j

    x_Width = x_Max - x_Min  REM 解析領域のx方向幅
    y_Width = y_Max - y_Min  REM 解析領域のy方向幅

    REM 0判定距離をセット
    If x_Width < y_Width Then
      Eps = y_Width * Eps0
    Else
      Eps = x_Width * Eps0
    End If

 

 このEpsの与え方は、解析領域の最大幅をスケール規模とするやり方です。これでは駄目な場合ももちろんあり得ますが、とりあえず(^^;)。今回は[メモ欄][宣言部][入力部]のみあげます。[計算部]は、引数並びから検討し直す必要があるとわかったので。

 

 ところでもう一個忘れてました。式(1)を見ると、定数πも必要です。これも[メモ欄][宣言部]に追加しておきます。

※ 後で調べてみたら、MAX関数とMIN関数と定数PIが組み込みになってました。まぁ~、好きにして下さい(^^;)


REM [メモ欄]

REM B_Count:境界要素数
REM EN_B(k,2):境界要素-節点対応表,k=1~B_Count
REM (k,1):要素kの始端の節点No.,(k,2):終端の節点No.

REM R_Count:領域要素数
REM EN_R(k,4):領域要素-節点対応表,k=1~R_Count
REM (k,j),j=1~4:jに従って要素kの要素節点No.を左回りに格納

REM BN_Count:境界節点数
REM BN(j,2):境界節点座標,j=1~BN_Count,(j,1):節点jのx座標,(j,2):y座標
REM RN_Count:領域節点数
REM RN(j,2):領域節点座標,j=1~RN_Count,(j,1):節点jのx座標,(j,2):y座標

REM BN_BP(j,2):節点の境界条件の有無,j=1~BN_Count
REM (j,1):節点jのψjが既知なら1、そうでないなら0
REM (j,2):節点jのqjが既知なら1、そうでないなら0

REM BN_BV(j,2):境界条件の値,j=1~BN_Count
REM (j,1):節点jのψjに与える値
REM (j,2):節点jのqjに与える値

REM z(i):未知ベクトル,i=1~2*BN_Count
REM BN_Z(j,2):節点-未知数対応表,j=1~BN_Count
REM (j,1):ψjのzの中の位置,(j,2):qjのzの中の位置
REM 節点j→ψj→z(j),節点j→qj→z(BN_Count+j)

REM A(i, j):全体係数マトリックス、i,j=1~2*BN_Count

REM Eps:0判定
REM Eps0:0オーダー

REM Pai:π


REM [宣言部]
REM 配列寸法の上限は1000

REM B_Count:境界要素数
REM R_Count:領域要素数
Dim EN_B(1000,2)  REM 境界要素-節点対応表,k=1~B_Count
Dim EN_R(1000,4)  REM 領域要素-節点対応表,k=1~R_Count

REM BN_Count:境界節点数
REM RN_Count:領域節点数
Dim BN(1000,2)  REM 境界節点座標,j=1~BN_Count
Dim RN(1000,2)  REM 領域節点座標,j=1~RN_Count

Dim BN_BP(1000,2) REM 節点の境界条件の有無,j=1~BN_Count
Dim BN_BV(1000,2) REM 境界条件の値,j=1~BN_Count

Dim BN_Z(1000,2) REM 節点-未知数対応表,j=1~BN_Count

Dim z(1000)  REM 未知ベクトル,i=1~2*BN_Count
Dim A(1000, 1000)  REM 全体係数マトリックス、i,j=1~2*BN_Count

REM Eps:0判定
Eps0 = 10^(-6)  REM 0オーダー

Pai = Acos(-1)  REM π:Fortranになじんだ人なら、懐かしくて泣いちゃう式(^^)


REM [入力部]

REM境界要素数B_Countと領域要素数R_Countは、事前に与えられると仮定する。
REM 境界要素と領域要素データは、要素ごとに節点番号が与えられると仮定する。
REM 与えられる節点番号は、左回りに順序付けられたものが与えられると仮定する。

REM 境界要素を要素番号kで読み込むLoop
REM  Loop内で要素ごとの節点番号を読み、境界要素-節点対応表EN_B(k,2)をセット

  For k = 1 to B_Count
    EN_B(k,1) = ?
    EN_B(k,2) = ?
  Next k

REM 領域要素を要素番号kで読み込むLoop
REM  Loop内で要素ごとの節点番号を読み、領域要素-節点対応表EN_R(k,4)をセット

  For k = 1 to R_Count
    For j = 1 to 4
      EN_R(k,j) = ?
    Next j
  Next k

REM境界節点数BN_Countと領域節点数RN_Countは、事前に与えられると仮定する。

REM 境界節点を節点番号jで読み込むLoop
REM   Loop内で境界節点座標を読み、BN(j,2)をセット

  For j = 1 to BN_Count
    BN(j,1) = ?
    BN(j,2) = ?
  Next j

REM 領域節点を節点番号jで読み込むLoop
REM   Loop内で領域節点座標を読み、RN(j,2)をセット

  For j = 1 to RN_Count
    RN(j,1) = ?
    RN(j,2) = ?
  Next j


REM 境界条件の場所を節点番号jで読み込むLoop
REM  Loop内で節点jの境界条件の有無を読み、BN_BP(j,2)をセット

  For j = 1 to BN_Count
    BN_BP(j,1) = ?
    BN_BP(j,2) = ?
  Next j

REM 境界条件の値を節点番号jで読み込むLoop
REM  Loop内で節点jの境界条件の値を読み、BN_BV(j,2)をセット

  For j = 1 to BN_Count
    BN_BV(j,1) = ?
    BN_BV(j,2) = ?
  Next j

REM 節点番号jのLoopで、節点-未知数対応表BN_Z(j,2)をセット
REM (j,1):ψjのzの中の位置,(j,2):qjのzの中の位置
REM 対応は、節点j→ψj→z(j),節点j→qj→z(BN_Count+j)

  For j = 1 to BN_Count
    BN_Z(j,1) = j
    BN_Z(j,2) = BN_Count + j
  Next j


REM 全ての入力完了 **********************************************
REM 解析領域のx方向幅とy方向幅を取得
  x_Min = 10^10    REM x座標最小値初期化
  x_Max = -10^10   REM x座標最大値初期化
  y_Min = 10^10    REM y座標最小値初期化
  y_Max = -10^10   REM y座標最大値初期化

  For j = 1 to BN_Count
    h = BN(j, 1)  REM j節点のx座標

    If h < x_Min Then
      x_Min = h
    End If

    If x_Max < h Then
      x_Max = h
    End If


    h = BN(j, 2)  REM j節点のy座標

    If h < y_Min Then
      y_Min = h
    End If

    If y_Max < h Then
      y_Max = h
    End If

  Next j

  x_Width = x_Max - x_Min  REM 解析領域のx方向幅
  y_Width = y_Max - y_Min  REM 解析領域のy方向幅

REM 0判定距離をセット
  If x_Width < y_Width Then
    Eps = y_Width * Eps0
  Else
    Eps = x_Width * Eps0
  End If

nice!(0)  コメント(0) 

有限要素法の選点法で、2階常微分方程式の境界値問題を解く [ネコ騙し数学]

有限要素法の選点法で、2階常微分方程式の境界値問題を解く

 

問題1 次の微分方程式の解を求めよ。

  fem1-001.png

【解】

この微分方程式の特性方程式は

  

したがって、この微分方程式の基本解はだから、一般解は

  

である。

境界条件より

  fem1-010.png

これを解くと、

  fem1-002.png

したがって、

  fem1-09.png

(解答終)

 

これを差分で解くのは芸がないということで、有限要素法選点法を用いてこの微分方程式を解くことにする。

それに先立ち、問題の微分方程式を

  

と置き換え、問題の微分方程式をzの微分方程式に書き換えると、

  fem1-003.png

になる。

こんとき、この微分方程式の境界条件は

  

となる。

そこで、この境界条件を満たす関数

  

を近似解と仮定するケロ。

そこで、(2)のzuに置き換えた

  fem1-004.png

を残差と呼ぶことにする。

もしuが(2)の解であれば、残差は0。しかし、(3)は(2)の解ではないので、R≠0である。

ここで、

  

となるように重み関数w(x)を調整し、(5)の積分の値が0になるようにする。

そして、(5)を満足する(3)を微分方程式(2)の近似解としようじゃないか。

 

この重み関数にδ関数

  

なるナゾの関数を採用し、(5)の積分を行うと

  

 

になる。

 

さてさて、

  

したがって、

  

x₀=1/2とすると、(6)より

  fem1-005.png

senten-1.pngだから、

 

で、

  

だから

  

に違いない。

 

赤線が厳密解の

  

青線が

  

この微分方程式に関する限り、厳密解との差は殆どないという驚きの結果が得られる。

 

1/3を選点x₀に選ぶと、

  fem1-006.png

したがって、

  

よって、

  

 

senten1-002.png

 

この微分方程式の場合、1/2にとったほうが1/3にとったほうが精度がよいことが分かる。

 

このように微分方程式の解を求める方法を選点法という。

 

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

  

【解】

同次方程式

  

の特性方程式は、

  

したがって、同次方程式の一般解は

  

特殊解y₀

  fem1-007.png

したがって、この微分方程式の一般解は

  

である。

境界条件x=0のときy=0より

  

x=1のときy=0より

  

よって、

  

(解答終)

 

この微分方程式の場合も同様に

  

とおくと、

  

となる。

そこで、x=1/2にとると、

  fem1-008.png

したがって、

  

 

senten1-003.png

 

粗い近似だから合うほうがおかしいケロよ。そこで、精度をあげる方法を考えることにするにゃ。

 


nice!(0)  コメント(0) 

ホイン(Heun)法 [ネコ騙し数学]

ホイン(Heun)法

 

§1 オイラー法

 

微分方程式

の解は

である。

(2)式の右辺の積分を

で近似すると、(2)式は

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

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

になる。

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

である。

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

 

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

【解】

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

x₁=0.1y₁=1.05だから

x₂=0.2y₂=1.163だから

したがって、

x=0.1のとき、y=1.05

x=0.2のとき、y=1.1106

x=0.3のとき、y=1.1846

(解答終)

 

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

 

 

 

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

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

 

 

§2 ホイン法

まず、オイラー法

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

と修正する。

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

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

 

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

 

Heun-001.png 

 

Heun-graph-001.png 

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

 

#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(予測修正子法ともいう)のほうが精度よく計算できていることがわかるだろう。


タグ:数値解析
nice!(0)  コメント(0) 

数値計算の誤差:打ち切り誤差と丸め誤差 [ネコ騙し数学]

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

 

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

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

有限項で打ち切れば、

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

 

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

かりに

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

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

と一致する。

 

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

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

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

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

と評価されるであろう。

 

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

たとえば、10進数の0.1は2進数で書くと

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

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

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

 

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

例えば、

という計算を考える。

いま仮に10進数で8桁(浮動小数)しか精度を有さない計算機があるとする。

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

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

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

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

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

と精度が保たれる。

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

 

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

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

 

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

となってしまう。

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

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

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

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

 

 

 

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

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

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

 

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

 


タグ:数値解析
nice!(0)  コメント(0) 

[境界要素法プログラムを設計する(基本のGauss掃き出し法)] [ネコ騙し数学]

[境界要素法プログラムを設計する(基本のGauss掃き出し法)]

 

1.Gauss掃き出し法の標準構成

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


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


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

 Gauss_0の引数Aは、解くべき連立一次方程式の係数マトリックスを表す配列,zGauss_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 ThenIfブロックはおいといて、REM文の後のサブルーティン呼び出しに注目です。呼ばれるサブルーティンは、スケーリング(Scaling),前進消去(Foward Elimination),後退代入(Back Substitution)を行います。このうち前進消去と後退代入の意味がわからない人は、

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

 

2.スケーリング

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

  

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

  

 (2')より、

  

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

  

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

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

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

 一方、

  

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

  

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

  

 (1")より、

  

 

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

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

 

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

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

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

 

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

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

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

 

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

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

 

 

 

3.ピボット選択

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

 

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

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

 

 

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

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

 

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

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

 

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

 

 スケーリングとピボット選択をセットで使うと、概ね6桁くらいの有効数字を数値誤差の影響から救えます。10^61,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²さん)



nice!(0)  コメント(0) 

極方程式で与えられた曲線に囲まれた部分の面積 [ネコ騙し数学]

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

 

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

  

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

 

kyoku-(x-a)^2+y^2=a^2.png例えば、半径aa>0)で点(a,0)を中心とする円

  

は、

  

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

 

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

  

であり、

  

という関係がある。

 

 

定理

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

  

である。

【略証】

  

(証明終)

 

 

Cardiod.png問題1 次の曲線(カージオイド)に囲まれている部分の面積を求めよ。

  

【解】

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

したがって、

  kyoku-men-001.png

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

  

(解答終)

 

Lemniscate.png問題2 次の曲線(レムニスケート)に囲まれている部分の面積を求めよ。

    

【解】

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

  

とおくと、

  

は、

  

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

  

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

  

よって、

  

したがって、

  

(解答終)

 

 

Clover3.png問題3 次の曲線(3葉線)に囲まれる部分の面積を求めよ。

  

【解】

とおくと、

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

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

  

になるのは、

  

したがって、

  kyoku-men-002.png

よって、

  

(解答終)

 


nice!(0)  コメント(0) 

テーラー級数、マクローリン級数を使って [ネコ騙し数学]

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

 

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

 

問題1

  

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

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

【解】

  

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

  

または、剰余項を用いて

  

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

  

にすればよいので

  

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

よって、

  

で、誤差は

  

(解答終)

 

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

  

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

 

 

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

  

の近似値を求めよ。

【解】

  

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

  

したがって、

  

よって、

  

(解答終)

 

 

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

  

【解】

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

だから、

  

とおくと、

  

よって、

m-siki-001.png

(解答終)

 

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

なお、

  


nice!(0)  コメント(0) 

差分法の基礎 [ネコ騙し数学]

差分法の基礎

 

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


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

ここで、記号Oはランダウの(ビッグ)O

  

ある。

 

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

  

が成り立つ。

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

  

を得ることができる。

 

(3)式より、h≠0のとき、

  

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

  sabun-004.png

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

 

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

  

より、h≠0のとき、

  

となり、f''(x)

  

で連続だから

  

となるMが存在し、

   

となり、

  sabun-004.png

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

sabun-graph-001.pngf(x)xの1次関数の時、

  

なので、

  

が成立すので、(6)は1次の精度である。

 

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


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

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

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

このことを、

  

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

 

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

  

が得られる。

sabun-graph-002.png(8)式で与えられるf'(a)の近似式を中心差分と呼び、この近似式の誤差はhの2次のオーダーであることと予想できる。

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

  

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

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

  

となるので、

  

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

 

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

 

の近似式を求めるために、(5)と(5’)の辺々を加えると、

  

sabun-graph-003.pngしたがって、

  

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

 

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


nice!(0)  コメント(0) 
前の10件 | 次の10件 ネコ騙し数学 ブログトップ