スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。


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

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

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。