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

Scilabでラグランジェ補間 その1ではforループを使って多項式補間を行うスクリプトを書きました。しかしながら、Scilabではループ計算は遅いので,できるだけ使わないのがコツなので、今回は多項式生成コマンドであるpolyを利用して高速化を行いました。(LagrangePoly_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

// *** 多項式(ラグランジュ)補間 ***
Wn = poly(0,"x");
Pn = 0;
for i = 1:length(Xn) do
Wn = poly(Xn(find(Xn <> Xn(i))),"x");
Pn = Pn + Fn(i) * Wn / horner(Wn,Xn(i));
end

PN = horner(Pn, X);

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

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



Scilabでループは遅い


Scilabでラグランジェ補間 その1では任意の関数に対する多項式補間を行う方法を書きました。その際に、多項式を計算するスクリプトとしてScilabで学ぶわかりやすい数値計算法サポートページにて公開されているlagrange.sceを参考にして、以下の様なコードを書きました。

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));


しかしながら、このコードはforが3つも入れ子になっています。
通常のプログラミング言語ではあまり問題にならないのかもしれませんが、Scilabに関してはforを多用すると計算速度が遅くなってしまいます。(参考:ループ計算は遅いので,できるだけ使わないScilabでは速いとされる)

そこで本ブログでもScilabで繰り返し計算: ロジスティック写像Scilabの論理演算で条件分岐などループを減らす工夫を行っています。

今回は、多項式の生成コマンドであるpolyを積極的に使って高速化を行いました。

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


N+1 個のデータ点 fn(n=1, 2, …, N+1) を全て通るような多項式は、以下の様に得られます。

数値計算の常識より

\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は"数値計算"を行うソフトウエアなので、数式を数式の形のままで微分したり、数式同士で割り算をしたり、因数分解したりということは基本的にはやりません。これらの仕事は本来Maximaなどの数式処理ソフトウエアの役目です。

しかしながら、全くそういった仕事ができないかと言うとそういうわけでもなく、多項式に対してこれらの処理を行うことが出来るコマンドを持っています。

今回利用したpolyは、多項式の解から多項式の係数を計算するコマンドです。(参考: scilab 根から多項式を作る!polyscilab 多項式の係数だけを取り出す)

逆に多項式の根を計算するにはrootsを使います。これはScilabで二酸化炭素の状態方程式 その2で利用しました。

多項式を「計算してしまって」数値にするのにはhornerを利用します。

関連エントリ




参考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シンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

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