風と翼 Windandwing は、c++の数値計算から風と翼の計算へ を目的とした数値計算を行っています。

Addres http://www5.synapse.ne.jp/windandwing/
c++ 数値計算 2  連立一次方程式( Jordan、Seidel、LU、変形 Cholesky )、                       補間法( Spline、Lagrange、 Aitken、直交関数法 )

第4章 連立一次方程式補間法(SplineLagrangeAitken直交関数法)

概要
1、連立一次方程式:
  ・Gauss-Jordan法では、STL/vectorを使用し、算出した。
  ・Gauss-Seidel法での計算実施。
  ・LU 分解法による解算出。
  ・変形コレスキー法による係数が対称行列である連立一次方程式の解算出。

2、補間法:
  Spline法(STL/vector使用、Gauss-Jordan法を適用)、Lagrange法、Aitken法を計算。
  直交関数であるLegendre多項式の関数近似法で、補間を試みた。
  なお、これらの補間法は、珪素鋼板のBH特性の補に適用し、精度や容易さを調査している。


1、連立一次方程式

1、Gauss-Jordan適用の連立一次方程式(コード1)
  STL/vector を使用し、Gauss-Jordan の解を算出した。通常、この計算は、2列の配列を用いて計算されるが、
  一列に数値が格納される STL/vector を用い、面白い計算を行った。

  この一列の vectorで、解を求める事ができる理由は、記述が2列であっても、数値の格納されている実際のアド
  レスは一列に並んでいる事を考慮して、STLの関数を用い、計算を行ったからである。

  実質のコードは30行で、行数は低減している。6元、2元、100元について計算を行った。
  この計算は、Spline補間に用いる。また、c++ 数値計算 8 行列式の逆行列の算出例でも用いた。

  行数低減
  max_element でアドレスを得る事で、行数が 2 行低減出来る。
  int inum=max_element(ddd.begin(),ddd.end())-ddd.begin();

2、Gauss-Seidel適用の連立一次方程式(コード2)
  STL/vectorを使用し、Gauss-Seidelで解を算出した。4元について計算を行った。

  Seidel法は、x1の初期値を決め、その初期値でx2の初期値を求め、更にx3の初期値を順次求める方法である。
  この方法に一括処理に優れたSTL/transform,accumulate等を適用し、簡素化するにはひと工夫が必要である。
  この例では、部分的にSTLを利用した。

3、LU分解による連立一次方程式(コード 2_1)
  LU分解法により、連立一次方程式の解の算出計算を行った。
  多くの本にこの計算方法について記載されているので、ここでは、説明しない。
  残念ながら STL の使用は vector のみである。しかし、n 行 n 列の n の次数によらず、データを csv ファ
  イルで読み込み、算出するようにした。現在、STL の活用できる文はないか、コードの再検討中 2015/4/2
  以下の文は、STL が活用出来そうだとして、後日検討する事にし、再掲する。2015/4/7
    for(int i=0;i<n;i++) intermediate_a[(n+1)*i]=1.0;

4、対称行列に対する変形コレスキーによる連立一次方程式(コード 2_2)
  対称行列(要素 aij が aij=aji となる)の連立一次方程式の解算出計算である。上三角行列と対角行列の要素の
  値を明示的に出力出来るようにコードを書いた。残念ながら c++ の STL の使用は vector のみである。
  しかし、n 行 n 列の n の次数によらず、データを csv ファイルで読み込み、算出するようにした。
  現在、STL の活用できる文はないか、コードの再検討中 2015/4/2
  以下の文は、STL が活用出来そうだとして、後日検討する事にし、再掲する。2015/4/7
   for(int i=0;i<n;i++) u[(n+1)*i]=1.0;

5、使用したSTLと関数
  単項関数オブジェクト(unary_function)、a.push_back()、swap()、a.assign()、a.insert()、transform()、
  max_element()、s.clear()、setw()、setprecision()、2項関数オブジェクト(binary_function)、
  accumulate()、inner_product()。

6、参考図書
  Fortranによる数値計算プログラム、STL標準講座、c言語と数値計算法、C による数値計算法、数値計算入門

2、補間法

1、Spline補間(コード3)
  多数項近似に比して、Spline補間は、Y軸の値の急激な変化に対応出来る補間法である。この補間法には
  Gauss-Jordanの連立一次方程式の解法を適用した。例は、Y軸の値が脈動したケースである。
   "cによる理工学問題の解法"に比してSTL使用によりコードの行数は少ない。

2Lagrange補間(コード4)
  Lagrange補間は、容易に数行で記述される補間法である。この補間法に、STLの適用は1行のみである。この補
  間法とSpline補間法の計算結果の比較を行う予定である。

3、Lagrange補間 STL多利用形(コード5)
  上記Lagrange補間に、STLを多く利用した計算法である。上記よりコードの行数が数行多くなった。
  STLのreplace()、accumulate()、inner_product() を使用している。

4、Lagrange補間 方程式の求根形(コード6)
 1)Lagrange補間は、y=an・x^n+a(n-1)・x^(n-1)+a(n-2)・x^(n-2)・・・a1・x+a0とし、xにx0を代入し
  て、補間値y0を求める方法である。この通常の方法で、xとyを入れ替えると、yに関する以下の方程式となる。

   bn・y^n+b(n-1)・y^(n-1)+b(n-2)・y^(n-2)・・・b1・y+b0-x0=0

  この多項式の bn は、上記のLagrange補間STL多利用形で算出し、yの求根は、c++STLによる方程式の求根
  二分法類似形の手法で算出した。

 2)評価 
  ・計算結果は、計算対象によるが、常時繰り返し回数が少なくて、高い精度を得る事が出来た訳ではない。
   これは、Lagrange補間(多項式近似)が持つy値がxのn乗のnの増加と伴に振動的に変化する特徴から生じて
   いると考える。(Runge 現象)

  ・一見すると、労多くして功少ない計算法であるが、y側で直接y値を求める方法は、補間が困難な域に適用で
   きるのではと期待している。コードの行数が多くなった。

5、Aitkenの逐次補間(コード7)
  Aitkenの補間は、2標本点間の一次補間多項式、3標本点間の2次補間多項式、、と次数を上げて、n次補間多項
  式と、n-1次との補間値差が指定したeps以下となった時に打ち切る補間法である。全標本点ではLagrange補間
  となる。
  計算過程で、STLの異母弟であるvalarrayを使用すれば、簡単になる個所もあるが、STLの各ライブラリーとの
  接続の際の煩雑さを考えると、STLで計算した。

  本補間は、f(x)が振動する計算例では、収束しない場合がある。

  注意点:補間点との距離(絶対値)で、同じ値が隣合せた場合、エラーが生じるので、補間点を1.0000001倍移動
      し、計算するようにした。

6、直交関数 Legendre 多項式による補間
 1)不均等間隔あるいは不等間隔の離散データを直交関数である Legendre 多項式で級数展開(数値計算 5 直交
   
関数 多項式級数展開)出来る。
   この事は、離散データを Legendre 多項式 Pn を用いて、補間出来る事を意味している。従い、以下が成立
  する。
     
                                  6-1式

 2)Legendre 多項式による補間計算を行うために、以下の計算手 順を考えた。

  (1) [ a,b ] 区間のデータの xk を [ -1,1 ] 区間 xj に変換する。これは、数値計算 6 数値積分とは、
    逆の変換である。

      
                                 6-2式

  (2) 位置 xj 毎に、適切な次数の Legendre 多項式を選定し xj を代入する。

  (3) 不等間隔用の適切な積分式に 各 xj ( j=0,1,...) および j の離散データを代入し、an を算出する。

  (4) 得た an で、6-1式により補間する。

  この補間の結果とコードは、BH曲線の補間法比較に示す。結論から言えば、BH曲線に対する適切な次数の
   Legendre 多項式、適切な積分式を得る事は、多くの計算を必要とした。


7、使用したSTLと関数
  unary_function、a.push_back()、swap()、a.assign()、a.insert()、transform()、*max_element()、
  setw()、setprecision()、u.resize()、adjacent_difference()、search_n()、greater<double>()。
  inner_production()、v.push_back()、x.resize()、v.assign()、replace()、accumulate()
  find_if()、min_elemnt()、fabs()、out.precision()、break、exit(1)、v.push_back()、x.resize()、
  transform()、sort()、break、less<double>()、findf()、fabs()、adjacent_difference()、count_if()、
  bind2nd(equal_to<double>()

8、参考図書
  cによる理工学問題の解法、STL標準講座、c言語と数値計算法

Information

     
2012年
この章の全体履歴
2012年1月5日 再公開 
2012年
この章の変更追加履歴
8)LU分解・変形コレスキー法の再掲 2015/4/7
7)変形コレスキー法の説明文を修正 2015/4/2
6)連立方程式の項に変形コレスキー法を追記 2015/4/1
5)連立方程式の項に LU 分解の計算を追記 2015/4/1
4)連立方程式の項に逆行列へのリンク挿入 2014/5/4
3)EVERGREENへの道等の写真の追加 2012/7/11
2)オペラ座・星・大地の写真の追加 2012/4/17
1)リンク先の誤記修正 2012/4/7

ナビゲーション


遠くに、みなとみらい

EVER GREEN への道

瀬戸大橋と船

ガーナ タマレ 大地の使いの家

 はるか彼方のオペラ
怪人ルパンとシャネルの女性のどちらもパリ
 遠いパリ 夜のオペラ座

 横浜みなとみらい 汽車道

環状2号線 シャネルの女性