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カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプPICCPA強磁性常微分方程式モンテカルロ解析odeトランジスタ状態密度インターフェーススイッチング回路ecaljPDS5022DOS定電流半導体シェルスクリプト乱数レベルシフトHP6632Aブレッドボード分散関係温度解析R6452Aトランジスタ技術I2C可変抵抗反強磁性セミナー数値積分確率論偏微分方程式バンド構造非線形方程式ソルババンドギャップ絶縁熱設計シュミットトリガLEDA/Dコンバータ三端子レギュレータLM358ISO-I2CGW近似カオスフォトカプラマフィンティン半径TL431数値微分PC817Cアナログスイッチ直流動作点解析発振回路USBサーボカレントミラー74HC4053パラメトリック解析LDAbzqltyチョッパアンプ量子力学FFT2ちゃんねるアセンブラBSch開発環境電子負荷ブラべ格子イジング模型補間基本並進ベクトル標準ロジック単振り子キュリー温度繰り返しMaxima状態方程式失敗談相対論スピン軌道相互作用FETランダムウォーク熱伝導六方最密充填構造コバルトewidthTLP621GGAQSGW不規則合金位相図抵抗SMPcygwinラプラス方程式スレーターポーリング曲線gfortranスイッチト・キャパシタ詰め回路TLP552三角波格子比熱TLP521条件分岐LM555MCUNE555QNAPマントルテスタ過渡解析FXA-7020ZRダイヤモンドデータロガーガイガー管自動計測Writer509UPSシュレディンガー方程式ブラウン運動awk差し込みグラフ熱力学平均場近似仮想結晶近似VCAfsolve井戸型ポテンシャルVESTA起電力スーパーセルOpenMP第一原理計算ubuntu固有値問題L10構造OPA2277interp12SC1815fccウィグナーザイツ胞面心立方構造フィルタジバニャン方程式ヒストグラム確率論マテリアルデザインspecx.f等高線正規分布PGAフェルミ面非線型方程式ソルバ初期値固定スピンモーメントスワップ領域ルチル構造リジッドバンド模型edeltquantumESPRESSO岩塩構造BaOSIC二相共存ZnOウルツ鉱構造フォノンデバイ模型c/aノコギリ波全エネルギーFSMTeXgnuplotmultiplotハーフメタルCapSense半金属合金結晶磁気異方性Ubuntu文字列入出力TS-110TS-112疎行列Excel直流解析ヒストグラム円周率不規則局所モーメントトラックボールPC等価回路モデルパラメータ・モデルキーボードRealforce三次元マンデルブロ集合フラクタル化学反応重積分縮退日本語最小二乗法関数フィッティングGimpMAS830LHiLAPW熱拡散方程式両対数グラフナイキスト線図負帰還安定性陰解法Crank-Nicolson法P-10クーロン散乱境界条件連立一次方程式片対数グラフEAGLEPIC16F785LMC662トランスシンボルCK1026線種凡例MBEAACircuitグラフの分割軸ラベルifort

最新コメント
リンク

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