Scilabで非線形方程式ソルバ その3

Scilabで非線形方程式ソルバ その1その2では関数の解を非線形方程式ソルバで計算しました。
今回は、離散データに対しても似たような解を求めるために、関数の補完と組み合わせるという事を行いました。

001_201502150526487fb.png
Fig.1: 離散データに対する補完と非線形方程式ソルバ



離散データのfsolve


Scilabで非線形方程式ソルバ その1その2では関数 f(x) がゼロとなる x を探す方法を書きました。今回は f(x) が具体的な関数の形ではなく、離散データとして与えられている場合について書きます。

と言ってもやっていることは簡単で、離散データを内挿(補間)してから解を求めるだけです。

今回はScilabで非線形方程式ソルバ その1との比較のために対数関数を例にしますが、実際の離散データは実験データであったり、複雑な数値計算(例えば第一原理計算や回路シミュレーション)の結果得られる数値列を想定します。

離散データの補完にはScilabでデータの補間で紹介したinterp1を利用します。
interp1の補間方法のデフォルトは線形補完(linear)のようです。今回のエントリでは最近接(nearest)は使えません。

スクリプト


スクリプトはlogsolve2_sce.txtとなりました。

今回は、離散データ自体もスクリプト内部で用意しましたが、実験データや、ほかの数値データを読み込む場合にはExcelデータを Scilabで読みこむScilabで大容量のCSV(テキスト)ファイルを読み込む,Scilabで大容量のCSV(テキスト)ファイルを読み込む2などを参照してください。

clear;

// *** データの作成 ***
X = linspace(0.1,10,15);
Y = log(X);

// *** 解くべき関数の定義 ***
function y = func(x)
// y = log(x) - yp;
// y = interp1(X, Y, x, "spline") - yp; // スプライン補完
y = interp1(X, Y, x, "linear") - yp; // 線形補完
endfunction

// *** 非線形方程式ソルバ ***
yp = 1; // yp = f(x)
x0 = 1; // ソルバ―の初期値
// 非線形方程式を解く
xp = fsolve(x0, func)
// 誤差
err = abs(xp - exp(yp))

// *** グラフのプロット ***
plot(X, Y, '.-b');
plot(X, yp * ones(Y), '--k');
plot(xp, yp, 'or');


結果


所詮、離散的な数値データから間の値を予測するだけなので、Scilabで非線形方程式ソルバ その1と比較すると結果に誤差が出ます。

線形補完を行った時の近似解と、真の値からの誤差は以下のようになりました。
 xp  =    2.737908
err = 0.0196261


スプライン補完を行った場合は、以下の通りです。
 xp  =    2.7130249
err = 0.0052569


今回の例ではスプライン補完の方が良い結果が得られているようです。

関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



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

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 非線形方程式ソルバ 補間 fsolve interp1 

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 補間 

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 補間 

Scilabでデータの補間

Scilabで数値積分:地球深部の密度と圧力では表として与えられているデータの数値積分を行いました。今回は表として与えられているデータの、データとデータの間の値を推定する、データの補間(内挿)の方法について書きます。

001_201311132310181fb.png

Fig.1: PREMの表に載っている地球の外核の密度(青丸)とその補間値(赤線)



Fig.1は、地球の半径1221.5kmから3480kmに位置する「外核」の密度です。
青の丸で示してあるのが、地震波観測から得られた一次元地球内部モデルPREM (Dziewonski and Anderson, 1981)です。PREMの表には半径200kmおきに密度などのデータが記載されています。しかしながら、例えば半径2500kmでの密度の値が知りたいとしても値が書いてありません。
こういうときには、どうすればよいでしょうか?いまの例だと、例えば、2400kmでの値と2600kmでの値が載っているので、平均を取るということも出来ます。

このように2点の間の点の値を2点間を結ぶ直線上の値として計算する方法を線形補間といいます。Scilabで数値積分:地球深部の密度と圧力で積分を計算するときに行ったのも実を言うとこれでした。

赤の実線でプロットされているのがPREMのデータからスプライン曲線で補間(内挿)をおこなったものです。
補間にはinterp1を利用しました。

[yp]=interp1(x, y, xp [, method,[extrapolation]])


x,yが保管される元となるデータでxpが補間したい値です。
xpは単一の数値でも、ベクトルでも大丈夫です。今回の例の様にxpにベクトルを入力するとypもベクトルで返ってきます。
今回の例ではスプライン補間をするため、methodの部分には'spline'を指定します。

補間の種類


Scilabの補間の方法には'spline','linear','nearest'の3種類が指定できます。
これらの違いを示すために3種類のプロットを行ったのがFig.2です。

002_20131113231018937.png

Fig.2: 3種類の補間の違い。スプライン補間(赤)、線形補間(青)、最近接データのプロット(緑)


赤のスプライン補間と青の線形補間はほとんど同じ値を示していますが、よく見るとスプライン補間の曲線のほうが上にt凸な形をしています。

少ないデータ点からグラフを滑らかに書くのには便利ですが、補間曲線の種類によって値が変わってしまうので値自体に過信は禁物です。

clear;
format('e',12);

// *** PREMのテーブルを読み出し ***
X = fscanfMat('PREM.txt');
Radius = 1E3 * X(:,1); // 半径 (m)
Vp = X(:,2); // P波速度 (m/s)
Vs = X(:,3); // S波速度 (m/s)
RHO = 1E3 * X(:,4); // 密度 (g/m^3)
Ks = 1E12 * X(:,5); // 断熱体積弾性率(Pa)
Mu = 1E12 * X(:,6); // 剛性率(Pa)
Nu = 1E12 * X(:,7); //
Pressure = 1E12 * X(:,8); // 圧力(Pa)
Gravity = X(:,9); // 重力加速度 (m/s^2)

// 半径と密度の外核の部分だけ取り出し
OCR = Radius(9:21); // 外核の半径 Outer Core Radius
OCD = RHO(9:21); // 外核の密度 Outer Core Density

// *** 内挿 ***
OCRp = linspace(OCR(1),OCR($));
OCDp = interp1(OCR,OCD,OCRp,'spline');

// *** プロット ***
// 外核全領域のプロット
scf(0);
plot(1E-3 * OCRp, 1E-3 * OCDp,'-r');
plot(1E-3 * Radius(9:21),1E-3 * RHO(9:21),'ob');
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");

// 2300-2700kmの領域のプロット
scf(1);
//xsetech([0,0,0.95,0.95]);
// スプライン曲線による補間
plot(1E-3 * OCRp, 1E-3 * interp1(OCR,OCD,OCRp,'spline'),'-r');
// 線形補間
plot(1E-3 * OCRp, 1E-3 * interp1(OCR,OCD,OCRp,'linear'),'--b');
// 最近接データの値による補間
plot(1E-3 * OCRp, 1E-3 * interp1(OCR,OCD,OCRp,'nearest'),'-.g');
plot(1E-3 * Radius(9:21),1E-3 * RHO(9:21),'ob');
zoom_rect([2300,11050,2700,11350]);
legend(['spline','linear','nearest']);
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



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

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 補間 

FC2カウンター
カテゴリ
ユーザータグ

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

最新コメント
リンク

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