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では、(解析的に解を得ることも可能な)非常に簡単な関数についてソルバ―を用いて解を求めました。

しかし前回の関数は、得られる解が一つだけでした。
では、解が二つ以上得られる場合、例えば二次関数の解を求める場合はどうなるのでしょうか?

002_20150212083236a8a.png

Fig.1: y = x2 と y = 1 の交点のx座標を計算した結果(上)。非線形方程式ソルバに与えた初期値x0とその結果として得られる解の関係(下)。


結論から言うと、非線形方程式ソルバでは解が一つしか得られず、その値はパラメータとして与える初期値に依存します。得られる解は、与えた初期値に近い値となりますが、必ずしも最も近い値になるというわけではないようなので、初期値の選び方は慎重に行う必要があります。


初期値の重要性


Scilabで非線形方程式ソルバ その1では、単純な対数関数の解を非線形方程式ソルバ―を用いて計算し、解析解と一致することを確認しました。

今回は x2 - 1 = 0 の解が、初期値 x0 の値によってどのように変化するかを確認します。

今回の記事で着目する二次関数と対数関数の違いは、得られる解が一つであるか二つであるかという事です。
Scilabの非線形方程式ソルバーは、数値的に解の近似値を計算するため、実際には複数の解が存在する場合であっても解が一つしか出てきません。
多くの場合、物理的に意味のある解は一つだけだったりするのですが、上手く重要な解を得るためには、初期値をどのように設定するかが重要になります。

x2 - 1 = 0 の解は当然ながら x=1 と x=-1 の二つなのですが、ソルバ―に与える初期値に応じてどちらの解が得られるかを確認するのが今回のエントリの趣旨です。

なお、多項式のすべての解をすべて一気に計算するための命令もScilabには用意されています。(参考: roots)

スクリプト


今回計算するSciabスクリプトはparabolicsolve_sce.txtです。

clear;

// *** データの作成 ***
X = linspace(-2,2);
Y = X .^ 2;

// *** 解くべき関数の定義 ***
function y = func(x)
y = x .^ 2 - yp;
endfunction

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

// *** 非線形方程式の初期値依存性 ***
X0 = X;
XP = fsolve(X0, func);

// *** グラフのプロット ***
subplot(2,1,1);
plot(X, Y, '-b');
plot(X, yp * ones(Y), '--k');
plot(xp, yp, 'or');
xlabel("x")
ylabel("y")

subplot(2,1,2);
plot(X0, XP, '-or')
xlabel("x0")
ylabel("xp")


結果


結果は、冒頭に紹介したFig.1です。
注目すべきは下のパネルの x0 = 0 近傍の挙動です。
必ずしも最も近い解が得られるわけではなく、初期値の微妙な変化で得られる解が振動していることがわかります。

なお、とにかくたくさん解がありそうな場合は、グラフを描いてみるのが良いです。Scilabでジバニャン方程式は f(x,y) = 0 となる x と y の組み合わせをすべてプロットしたものです。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 非線形方程式ソルバ fsolve 初期値 

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

最新コメント
リンク

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