Scilabでラグランジェ補間 その1

Scilabでデータの補間ではScilabの標準のコマンドであるinterp1を用いて線形補間やスプライン補間などを行いました。
今回は、補間の考え方の勉強をするために、(実際にはあまり使われない)ラグランジェ補間(多項式補間)をやってみました。

001_20140629231048c59.png

Fig.1: 関数f(x)=2exp(x-1)-1 (青実線)とその3点(赤四角)を通るように計算した多項式補間(黒破線)、及び関数と補間値の誤差(下パネル)。



補間(内挿)とは


数値計算の常識には、内挿に関して以下の様に書かれています。(ただしnの番号付けだけはScilabでスクリプトを書くことを考慮し0からではなく1からに変更しました。これ以降も同様です。)

関数f(x)の値が, 変数xのとびとびの異なる値 x1, x2, …, xN+1に対して与えられているとき, すなわち

fn ≡ f(xn) (n = 1,2, …, N+1)

のみが与えられているとき, xn以外のxに対するf(x)の値を"推測"しようというのが広い意味での"補間"(interpolation; "内挿"ともいう)である.


補間はいろいろな場面で利用されます。例えばScilabでデータの補間AkaiKKRでコバルトの格子定数などです。

今回も例によって数値計算の常識を教科書に補間について勉強します。数値計算の常識の補間に関する大雑把な内容は、以下の通りです。

  1. 最初に考え付くのは多項式(ラグランジェ)補間。しかし、上手く行かない場合も多い。
  2. 次に思いつくのが線形補間。でも折れ曲がりが気になる。
  3. 結局良く使われているのがスプライン補間。


2と3の線形補間とスプライン補間に関してはScilabに標準の命令がありScilabでデータの補間で紹介しています。
そこで今回は(あまり使われない)ラグランジェ補間について書きます。

多項式補間(ラグランジェ補間)


多項式補間の基本的な考え方は N+1 個のデータ点 fn(n=1, 2, …, N+1) を全て通るようなN次多項式は必ず一意に決まるはずなので、それを近似式としましょうと言うものです。

数値計算の常識より

\omega (x) = (x - x_1)(x - x_2) \cdots (x - x_{N+1})

\omega_n (x) = \frac{\omega (x)}{x - x_n}

とおくと、求める多項式 PNは以下の様に求められます。

P_N (x) = \sum_{n=1}^{N+1} \frac{\omega_n (x)}{\omega_n (x_n)} f_n

Scilabスクリプト


Scilabで学ぶわかりやすい数値計算法サポートページにて公開されているlagrange.sceを流用して以下の様なスクリプトを作成しました。(Lagrange_sce.txt)
clear;

// *** 近似する元の関数 ***
function y = f(x)
// y = 2 .* exp(x - 1) - 1
y = 2 ./ (1 + 9 .* x .^ 2) - 1
endfunction

// *** データ点 ***
N = 10; // データ数
Xn = linspace(-1, 1, N + 1); // xのデータ点
Fn = f(Xn); // yのデータ点
X = linspace(-1, 1, 100); // プロット用のx

// *** 多項式(ラグランジュ)補間 ***
Pn = zeros(X);
for m = 1:length(X)
for i = 1:(N + 1)
Z(i) = 1;
for j = 1:(N + 1)
if i <> j
Z(i) = Z(i) * (X(m) - Xn(j)) / (Xn(i) - Xn(j));
end
end
Pn(m) = Pn(m) + Fn(i) * Z(i);
end
end

// *** グラフのプロット ***
subplot(2,1,1);
plot(X, f(X));
plot(X, Pn, '--k');
plot(Xn, Fn, 'sr');

// *** 誤差のプロット ***
subplot(2,1,2);
plot(X, Pn - f(X));


結果


近似する関数は数値計算の常識に倣って以下の2種類を行いました。

f(x)=2 \exp(x-1) - 1

f(x)=\frac{2}{1+9x^2}-1

定義域はともに -1 ≦ x ≦ 1 です。
以下に示すように、前者は上手く補間できるのですが、後者はイマイチです。

002_20140629231048cb2.png

Fig.2: 関数f(x)=2exp(x-1)-1 (青実線)とその7点(赤四角)を通るように計算した多項式補間(黒破線)、及び関数と補間値の誤差(下パネル)。Fig.1と比較して、通るべき点の数を増やすことによって近似の誤差が少なくなっていることが読み取れる。


003_20140629231048721.png

Fig.3: 関数 2/(1+9x2)-1 (青実線)とその11点(赤四角)を通るように計算した多項式補間(黒破線)、及び関数と補間値の誤差(下パネル)。通るべき点の数を増やしても近似の誤差が小さくならない問題点がある。


Fig.1-2は補間が上手く行っている例です。しかしながら、Fig.3の様に補間するもとの関数によっては誤差が大きくなってしまうことがあるのがラグランジェ補間の問題点です。

スプライン補間


もう少し優等生なスプライン補間を用いた結果も比較のため載せておきます。(spline_sce.txt)

004_20140629233541671.png

Fig.4: スプライン補間。図の読み方はFig.3と同じ。Fig.3と比較して誤差が少なくなっていることが分かる。


Fig.3とFig.4を比較するとスプライン補間のほうが良い結果を示していることがわかります。
ただ、「スプライン補間さえ使ってればいつでもオッケー」と言うわけではないと言うことも覚えておかなければいけません。

数値計算の常識曰く

どだい, f(xn) (n = 1, 2, …, N+1) から f(x) (x≠xn)を知ろうというのは, 関数fがよっぽど良い性質を持っていなければ"原理的に"無理であるのは明らかであろう. (x = x1, …, xN+1 で一致するけれど他のxでは値が異なる関数はいくらでもありうる.)


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。


tag: Scilab 補間 

comment

Secret

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj常微分方程式モンテカルロ解析状態密度トランジスタodeDOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト分散関係レベルシフト乱数HP6632AR6452A可変抵抗トランジスタ技術温度解析ブレッドボードI2C反強磁性確率論数値積分セミナーバンドギャップバンド構造偏微分方程式非線形方程式ソルバ熱設計絶縁ISO-I2Cカオス三端子レギュレータLM358GW近似マフィンティン半径A/DコンバータフォトカプラシュミットトリガLEDPC817C発振回路数値微分直流動作点解析サーボカレントミラーTL431アナログスイッチUSB74HC4053bzqltyVESTA補間電子負荷アセンブライジング模型BSch量子力学単振り子2ちゃんねるチョッパアンプLDA開発環境基本並進ベクトルFFT標準ロジックブラべ格子パラメトリック解析抵抗SMPMaxima失敗談ラプラス方程式繰り返し位相図スイッチト・キャパシタ熱伝導状態方程式キュリー温度gfortranコバルトTLP621不規則合金Quantum_ESPRESSO六方最密充填構造ランダムウォーク相対論ewidthスピン軌道相互作用FETQSGWVCAcygwinスレーターポーリング曲線GGA仮想結晶近似PWscfシュレディンガー方程式LM555ハーフメタル固有値問題NE555最小値ガイガー管QNAPUPS自動計測ダイヤモンドマントルTLP552格子比熱最適化MCU井戸型ポテンシャル最大値xcrysdenCIF条件分岐詰め回路フェルミ面差し込みグラフスーパーセルfsolveブラウン運動awk過渡解析起電力三角波第一原理計算FXA-7020ZRWriter509Ubuntuテスタ熱力学データロガーTLP521OpenMPubuntu平均場近似MAS830LトランスCK1026PIC16F785PGA2SC1815EAGLEノコギリ波負帰還安定性ナイキスト線図MBEOPA2277P-10フィルタCapSenseAACircuitLMC662文字列固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型電荷密度重積分SIC二相共存磁気モーメント不純物問題PWgui擬ポテンシャルゼーベック係数ZnOウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデルcif2cell入出力陰解法熱拡散方程式HiLAPW両対数グラフCrank-Nicolson法連立一次方程式specx.fifort境界条件片対数グラフグラフの分割円周率ヒストグラム不規則局所モーメントGimpシンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ