スポンサーサイト

上記の広告は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カウンター
カテゴリ
ユーザータグ

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ブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolveCIFxcrysden最大値最小値ubuntu最適化平均場近似OpenMPシュレディンガー方程式固有値問題井戸型ポテンシャル2SC1815TeX結晶磁気異方性OPA2277非線型方程式ソルバフラクタルFSM固定スピンモーメントc/agnuplotPGA全エネルギーfccマンデルブロ集合縮退正規分布キーボード初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存ZnOウルツ鉱構造BaOSIC重積分磁気モーメント電荷密度化学反応クーロン散乱岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルRealforceフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線凡例軸ラベル線種シンボルトラックボールグラフの分割MAS830LPIC16F785トランス入出力CK1026PC直流解析パラメータ・モデル等価回路モデル不規則局所モーメント関数フィッティング日本語ヒストグラムTS-112ExcelGimp円周率TS-110LMC662片対数グラフ三次元specx.fifortUbuntu文字列疎行列不純物問題ジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式AACircuit熱拡散方程式HiLAPW両対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE最小二乗法

最新コメント
リンク

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