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

最新コメント
リンク

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