Scilabで最急降下法 その2

Scilabで最急降下法 その1では1変数の関数 f(x) に対して最急降下法を用いて最小値を求めました。今回は2変数の関数 f(x,y) について最小値を求めます。

001_20170509143114675.png
Fig.1: 最急降下法による最小値探査



具体的には二次元のガウス関数に-1を掛けた関数を f(x,y) とします。
\begin{equation}
f(x,y)=\frac{1}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}}\exp\left\{-\frac{1}{2(1-\rho^2)}\\
\times\left(\frac{(x-\mu_1)^2}{\sigma_1^2}-2\rho\frac{(x-\mu_1)(y-\mu_2)}{\sigma_1\sigma_2} +\frac{(y-\mu_2)^2}{\sigma_2^2}\right) \right\}
\end{equation}
ここで
\begin{equation}
\rho=\frac{\sigma_{12}}{\sigma_1\sigma_2}
\end{equation}

最急降下法の値の更新は二次元の場合は以下のように行います。
\begin{equation}
\begin{pmatrix}
x^{(k+1)}\\
y^{(k+1)}
\end{pmatrix}
=
\begin{pmatrix}
x^{(k)}\\
y^{(k)}
\end{pmatrix}
-\alpha
\begin{pmatrix}
\frac{\partial f(x^{(k)}, y^{(k)})}{\partial x}\\
\frac{\partial f(x^{(k)}, y^{(k)})}{\partial y}
\end{pmatrix}
\end{equation}
停止条件は ∂f/∂x < ε かつ ∂f/∂y < ε とすればよいと思います。

clear;

// *** 二次元ガウス分布 ***
mu1 = 3; mu2 = 1;
mu = [mu1; mu2];
sigma1 = sqrt(10); sigma2 = sqrt(10); sigma12 = 5;
SIGMA = [sigma1^2, sigma12; sigma12, sigma2^2];
rho = sigma12 / (sigma1 * sigma2);
function z = func(x, y)
z = -1 ./ (2 * %pi * sigma1 * sigma2 * sqrt(1 - rho)) ..
.* exp(-1 / (2 .* (1 - rho^2)) ..
.* ((x - mu1) .^ 2 ./ (sigma1 ^ 2) - 2 .* rho .* (x - mu1) .* (y - mu2) ./ (sigma1 * sigma2) + ((y - mu2) .^ 2) ./ (sigma2 ^ 2)))
endfunction


// *** 数値微分 ***
h = 1E-3;
function dzx = dfx(x, y)
dzx = (func((x+h), y) - func((x-h), y)) ./ (2 * h)
endfunction

function dzy = dfy(x, y)
dzy = (func(x, (y+h)) - func(x, (y-h))) ./ (2 * h)
endfunction

// *** グラフのプロット ***
x = linspace(-10,10);
y = linspace(-10,10);
[X,Y] = ndgrid(x,y);
Z = func(X,Y);
// 色の設定
//xset("colormap",jetcolormap(64))
// xset("colormap",graycolormap(64))
// Sgrayplot(x,y,Z)
xset("fpf"," "); // 等高線に値を表示しない
contour2d(x,y,Z,10);
xmin = min(X); xmax = max(X); ymin = min(Y); ymax = max(Y);
zoom_rect([xmin,ymin,xmax,ymax]);




// *** 最小値の計算 ***
// 停止条件
err = 1E-5;
a = 200;
// 初期値
x = -4; y = 4;
z = func(x, y);
dx = dfx(x, y);
dy = dfy(x, y);
plot(x, y, "xk");
// 最小値の計算
while ((abs(dx) > err) | (abs(dy) > err))
x = x - a * dx;
y = y - a * dy;
dx = dfx(x, y);
dy = dfy(x, y);
plot(x, y, ".r");
end
// 計算結果
x, y
plot(x, y, "xk");


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 最適化 最小値 最大値 

Scilabで最急降下法 その1

Scilabで何らかの関数 f(x) の最小値(あるいは最大値)を計算することを考えます。関数の値を計算するのが簡単な場合は x の定義域全体で f(x) を計算した後 minmax を使うという方法もあります。しかしながら f(x) の計算にそれなりの時間がかかる場合や f(x, y) といったように引数がたくさんある場合は効率的ではないと思います。

そこで今回は最急降下法のアルゴリズムを利用して f(x) の最小値を求めるということをやってみます。

001_20170507020049254.png

Fig.1: 最急降下法での最小値探索。上が関数f(x)の値、下が微分値f'(x)



最小値を求める関数


さて、実際に最小値を求める関数 f(x) ですが、今回は単純にガウス関数に -1 を掛けたものにします。
\begin{equation}
f(x) = - \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2\sigma^2}\right)
\end{equation}
当然ながら、f(x)が最小になるのは x = μ のときです。

最急降下法


高校の数学で習ったとおり f(x) が最大値や最小値(や極値)をとるときその微分は f'(x) = df(x)/dx = 0 となります。最急降下法は、関数の微分を計算しその傾きが大きいほうへ f'(x)=0 となる x を探すアルゴリズムです。具体的には以下の手続きを繰り返します。
  1. x の初期値 x(0) を決める
  2. f'(x) < ε なら終了
  3. x(k+1) = x(k) - αf'(x(k))
  4. 2.に戻る

実際には α や ε を上手に決めておく必要があります。αは勾配の方向にどの程度進むかを決めるパラメータ(下記Scilabスクリプトではa)で、εは計算の終了条件を決めるパラメータ(下記Scilabスクリプトではerr)です。

Scilabスクリプト


clear;

// *** 一次元ガウス分布 ***
function y = func(x)
mu = 3;
sigma = 1;
y = -1 / sqrt(2*%pi*sigma^2) * exp(-1 * ((x - mu) .^ 2) ./ (2*sigma^2))
// y = cos(x)
endfunction

// *** 数値微分 ***
function y = dfunc(x)
h = 1E-4;
y = (func(x+h) - func(x-h)) ./ (2*h)
endfunction

// *** グラフのプロット ***
X = linspace(0,6);
Y = func(X);
dY = dfunc(X);
subplot(2,1,1);
plot(X, Y);
subplot(2,1,2);
plot(X, dY);

// *** 最小値の計算 ***
// 停止条件
err = 1E-3;
a = 0.5;
// 初期値
x = 1;
y = func(x);
dx = dfunc(x);
subplot(2,1,1);
plot(x, y, ".r");
subplot(2,1,2);
plot(x, dx, ".r");

// *** 最小値の計算 ***
while abs(dx) > err
x = x - a * dx;
dx = dfunc(x);
y = func(x);
subplot(2,1,1);
plot(x, y, ".r");
subplot(2,1,2);
plot(x, dx, ".r");
end

// *** 計算結果 ***
x
subplot(2,1,1);
plot(x, y, "xk");


参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 最適化 最小値 最大値 

Scilabで二重積分

Scilabを利用すると1変数の数値積分が簡単に計算できます。

\begin{equation}
\int_{x_0}^{x_1}f(x)\mathrm{d}x
\end{equation}

このブログでも数値積分タグにいくつかの例を見つけることができます。しかしながら、2変数の数値積分はこれまで行ってきませんでした。

\begin{equation}
\int\int f(x,y) \mathrm{d}x\mathrm{d}y
\end{equation}

Scilabには二重積分を計算することが可能な int2d が存在します。今回は高校数学の美しい物語で解析的に解かれている二重積分を数値的に計算してみます。


積分範囲が長方形領域の場合


積分範囲が長方形の領域の場合、すなわち以下のような式で表すことができる場合は、簡単に数値積分できます。

\begin{equation}
\int_{x_0}^{x_1}\int_{y_0}^{y_1}f(x,y)\mathrm{d}x\mathrm{d}y
\end{equation}

Scilabの int2d では長方形領域を2つの三角形のパッチワークとして与えます。
積分範囲を int2d に渡すために行列XとYを用意します。それぞれ2つの三角形の頂点のx座標とy座標を与えます。

\begin{equation}
X =
\begin{pmatrix}
x_{0} & x_{0} \\
x_{1} & x_{1} \\
x_{1} & x_{0} \\
\end{pmatrix},
Y =
\begin{pmatrix}
y_{0} & y_{0} \\
y_{0} & y_{1} \\
y_{1} & y_{1} \\
\end{pmatrix}
\end{equation}

001_20170423145553b2b.png
Fig.1: Scilabのint2dへの積分範囲の与え方


実際に以下の積分を計算して見ます。

\begin{equation}
\int_{0}^{\pi}\int_{0}^{R}x^4 \sin(y)\mathrm{d}x\mathrm{d}y
\end{equation}

clear;

r = 1;

// *** 積分する関数の定義 ***
function z = func(x,y)
z = (x .^ 4) * sin(y)
endfunction
// 積分範囲
x0 = 0; x1 = r;
y0 = 0; y1 = %pi;

// *** 二重積分 ***
X = [x0, x0;
x1, x1;
x1, x0];
Y = [y0, y0;
y0, y1;
y1, y1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = 2*(r^5)/5


数値化解と解析解が同じ値になることが確認できます。

積分範囲が三角形の組み合わせで表せる場合


積分範囲が長方形の場合は2つの三角形の組み合わせで表現されますが、より複雑な形状の場合も任意の個数の三角形の組み合わせで表現できるはずです。今回は逆に簡単になってしまいますが、1個の三角形で表現できる例を計算します。

\begin{equation}
\int \int_D xy^2 \mathrm{d}x\mathrm{y}
\end{equation}

jusekibun.png
Fig.2: 積分領域Dが三角形ひとつ分の例


積分領域が三角形ひとつ分なので、与える行列は3行1列になります。

\begin{equation}
X =
\begin{pmatrix}
x_{0} \\
x_{1} \\
x_{1} \\
\end{pmatrix},
Y =
\begin{pmatrix}
y_{1} \\
y_{0} \\
y_{1} \\
\end{pmatrix}
\end{equation}

この計算を行うScilabスクリプトは以下のようになります。

clear;

// *** 積分する関数の定義 ***
function z = func(x,y)
z = x .* (y .^ 2)
endfunction

// *** 二重積分 ***
X = [0;
1;
1];
Y = [1;
0;
1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = 3/20


このスクリプトも数値解と解析解が同じに値になることが分かります。

同様にしてN個の三角形の組み合わせで表現される積分範囲の場合3行N列の行列で指定することができます。

更に複雑な積分領域の場合


どんなに複雑な積分領域の形状であっても三角形のパッチワークで表現できるはずですが、現実的には大変です。そこでOctaveの精義―フリーの高機能数値計算ツールを使いこなすで紹介されている方法を試してみましたが、現状うまく行っていません。上手く行っていませんがとりあえず方法だけは紹介します。
具体的にはScilabの論理演算で条件分岐の考え方を使って積分領域外では値がゼロになるように被積分関数の定義を行います。

\begin{equation}
\int\int_D -\frac{1}{(2x + y + 1)^2}\mathrm{d}x\mathrm{d}y
\end{equation}

jusekibun2.png
Fig.3: 複雑な積分領域の例


clear;

// *** 積分する関数の定義 ***
function z = func(x,y)
region = y >= x .^ 2
z = - 1 ./ ((2 * x + y + 1) .^ 2) .* region
endfunction

// *** 二重積分 ***
X = [0;
1;
1];
Y = [0;
0;
1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = (1/3) * log(4) - 1/2


原理的にはこのスクリプトでよいはずですが、実際には正しく計算してくれません。Scilab 6.0ではエラーで停止します。Scilab 5.5.2ではそれっぽい値を返しますが、解析解の値とはかなりずれた値となっており、不正確です。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 数値積分 重積分 

Scilabで数値微分 その4

Scilabで数値微分 その1, その2で行った誤差の議論を、誤差の最大値と誤差の平均値について行いました。

こういう議論に意味があるかはわかりませんが。


Scilabで数値微分 その1, その2ではsin(x)のx=0.3πのときの微分値を数値計算するに際して、差分の刻み幅を変えた際に数値計算の誤差(丸め誤差、打切り誤差)がどのように影響するかについて調べました。
今回はxを1から2πの範囲で微分値を計算し、差分の刻み幅を変えた際に数値計算の誤差の平均値や最大値への影響を調べます。

001_20150405230943e9c.png
Fig.1: 数値微分と解析解の比較(丸は誤差の平均値、四角は誤差の最大値)

002_201504052309435a9.png
Fig.2: 数値微分の刻み幅を2倍にした時の変化の度合い(丸は誤差の平均値、四角は誤差の最大値)


Scilabのスクリプトはdiff3.sceです。Scilabで数値微分 その3で作成したdifferential.sciが同じディレクトリにある必要があります。

clear;

// *** 関数の定義を読み出し ***
exec('differential.sci',0);

// *** 計算の設定 ***
xmin = 0;
xmax = 2 * %pi;
x = linspace(xmin,xmax);
dy = cos(x); // 解析解
N = [0:1:50]'; // 刻み幅

// *** 数値微分 ***
for i = 1:length(N)
// 刻み幅
n = N(i);
dx = 1 / 2 ^ n;
// 前進差分
dyf1 = diff_f1(x, dx, sin);
dyf12 = diff_f1(x, 2 * dx, sin);
ERRmeanf1(i) = mean(abs(dy - dyf1));
ERRmaxf1(i) = max(abs(dy - dyf1));
DIFmeanf1(i) = mean(abs(dyf12 - dyf1));
DIFmaxf1(i) = max(abs(dyf12 - dyf1));
// 中心差分
dyf2 = diff_f2(x, dx, sin);
dyf22 = diff_f2(x, 2 * dx, sin);
ERRmeanf2(i) = mean(abs(dy - dyf2));
ERRmaxf2(i) = max(abs(dy - dyf2));
DIFmeanf2(i) = mean(abs(dyf22 - dyf2));
DIFmaxf2(i) = max(abs(dyf22 - dyf2));
// 前進差分に対するRomberg1段
dyf1r = diff_f1r(x, dx, sin);
dyf1r2 = diff_f1r(x, 2 * dx, sin);
ERRmeanf1r(i) = mean(abs(dy - dyf1r));
ERRmaxf1r(i) = max(abs(dy - dyf1r));
DIFmeanf1r(i) = mean(abs(dyf1r2 - dyf1r));
DIFmaxf1r(i) = max(abs(dyf1r2 - dyf1r));
// 中心差分に対するRomberg1段
dyf2r = diff_f2r(x, dx, sin);
dyf2r2 = diff_f2r(x, 2 * dx, sin);
ERRmeanf2r(i) = mean(abs(dy - dyf2r));
ERRmaxf2r(i) = max(abs(dy - dyf2r));
DIFmeanf2r(i) = mean(abs(dyf2r2 - dyf2r));
DIFmaxf2r(i) = max(abs(dyf2r2 - dyf2r));
end

// *** グラフのプロット ***
// *** 誤差のプロット ***
scf(0);
a = gca();
a.data_bounds = [min(N),1E-14; max(N),1];
a.log_flags = "nl";
// 誤差の平均値
plot(N, ERRmeanf1, '-or'); // 前進差分
plot(N, ERRmeanf2, '-om'); // 中心差分
plot(N, ERRmeanf1r, '-ob'); // 前進差分に対するRomberg1段
plot(N, ERRmeanf2r, '-og'); // 中心差分に対するRomberg1段
// 誤差の最大値
plot(N, ERRmaxf1, '-sr'); // 前進差分
plot(N, ERRmaxf2, '-sm'); // 中心差分
plot(N, ERRmaxf1r, '-sb'); // 前進差分に対するRomberg1段
plot(N, ERRmaxf2r, '-sg'); // 中心差分に対するRomberg1段
// *** 刻み幅を変えた際の値の変化 ***
scf(1);
a = gca();
a.data_bounds = [min(N),1E-14; max(N),1];
a.log_flags = "nl";
// 差の平均値
plot(N, DIFmeanf1, '-or'); // 前進差分
plot(N, DIFmeanf2, '-om'); // 中心差分
plot(N, DIFmeanf1r, '-ob'); // 前進差分に対するRomberg1段
plot(N, DIFmeanf2r, '-og'); // 中心差分に対するRomberg1段
// 差の最大値
plot(N, DIFmaxf1, '-sr'); // 前進差分
plot(N, DIFmaxf2, '-sm'); // 中心差分
plot(N, DIFmaxf1r, '-sb'); // 前進差分に対するRomberg1段
plot(N, DIFmaxf2r, '-sg'); // 中心差分に対するRomberg1段


関連エントリ




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 数値微分 

Scilabで数値微分 その3

Scilabで数値微分 その1で紹介した数値微分の近似式をScilabのfunctionにまとめました。

// *** 前進差分 ***
function df = diff_f1(x, dx, func)
df = (func(x + dx) - func(x)) ./ dx
endfunction

// *** 中心差分 ***
function df = diff_f2(x, dx, func)
df = (func(x + dx) - func(x - dx)) ./ (2 * dx)
endfunction

// *** 前進差分に対するRomberg1段 ***
function df = diff_f1r(x, dx, func)
df = 2 * (diff_f1(x, dx, func) - 0.5 * diff_f1(x, 2 * dx, func))
endfunction

// *** 中心差分に対するRomberg1段 ***
function df = diff_f2r(x, dx, func)
df = (4 / 3) * (diff_f2(x, dx, func) - 0.25 * diff_f2(x, 2 * dx, func))
endfunction



上記の内容をdifferential.sciのように別ファイルにまとめておけば、以下のように呼び出して使えます。

clear;

// *** 関数の定義を読み出し ***
exec('differential.sci',0);

// *** 計算範囲と刻み幅 ***
n = 10;
xmin = 0;
xmax = 2 * %pi;
dx = 1 / 2 ^ n;
x = linspace(0,xmax);

// *** 正弦波とその微分 ***
y = sin(x);
//dy = diff_f1(x, dx, sin); // 前進差分
//dy = diff_f2(x, dx, sin); // 中心差分
//dy = diff_f1r(x, dx, sin); // 前進差分に対するRomberg1段
dy = diff_f2r(x, dx, sin); // 中心差分に対するRomberg1段

// *** グラフのプロット ***
// グラフの描画
plot(x, y, '-b');
plot(x, dy,'--r');
// グラフの装飾
legend(["sin(x)","d(sin(x))/dx"],3);
zoom_rect([xmin,-1,xmax,1]);


数値微分なので、微分値が具体的な数値として求まります。
xには微分したいxの値を、dxには微小な値を入力しますが、具体的にどの程度の数値が必要なのかはScilabで数値微分 その2を参考に決めます(hの値がこのエントリにおけるdxに相当します)。

そもそも関数の解析的な積分はできない場合が多い反面、微分はできてしまうことが多いので、数値的に微分値を求める需要はあまり多くないのかもしれませんが。

関連エントリ




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 数値微分 

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

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性OPアンプPICCPA常微分方程式モンテカルロ解析ecaljodeトランジスタ状態密度インターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプトレベルシフト乱数HP6632AR6452AI2C可変抵抗分散関係トランジスタ技術ブレッドボード温度解析反強磁性確率論バンドギャップセミナー数値積分熱設計非線形方程式ソルババンド構造絶縁偏微分方程式ISO-I2CLM358フォトカプラ三端子レギュレータカオスLEDシュミットトリガGW近似A/Dコンバータ発振回路PC817C直流動作点解析USBマフィンティン半径数値微分アナログスイッチTL43174HC4053カレントミラーサーボ量子力学単振り子チョッパアンプ補間2ちゃんねる開発環境bzqltyFFT電子負荷LDAイジング模型BSch基本並進ベクトルブラべ格子パラメトリック解析標準ロジックアセンブラ繰り返し六方最密充填構造SMPコバルトewidthFET仮想結晶近似QSGW不規則合金VCAMaximaGGA熱伝導cygwinスレーターポーリング曲線キュリー温度スイッチト・キャパシタ失敗談ランダムウォークgfortran抵抗相対論位相図スピン軌道相互作用VESTA状態方程式TLP621ラプラス方程式TLP552条件分岐NE555LM555TLP521マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測QNAPUPSWriter509ダイヤモンドデータロガー格子比熱熱力学awkブラウン運動起電力スーパーセル差し込みグラフ第一原理計算フェルミ面fsolve最大値xcrysden最小値最適化ubuntu平均場近似OpenMP井戸型ポテンシャルシュレディンガー方程式固有値問題2SC1815結晶磁気異方性OPA2277非線型方程式ソルバTeXgnuplot固定スピンモーメントFSMPGAc/a全エネルギーfccフラクタルマンデルブロ集合正規分布縮退初期値interp1multiplotフィルタ面心立方構造ウィグナーザイツ胞L10構造半金属二相共存SICZnOウルツ鉱構造BaO重積分クーロン散乱磁気モーメント電荷密度化学反応CIF岩塩構造CapSenseノコギリ波デバイ模型ハーフメタルキーボードフォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt合金等高線線種凡例シンボルトラックボールPC軸ラベルグラフの分割トランス文字列CK1026MAS830L直流解析Excel不規則局所モーメントパラメータ・モデル入出力日本語最小二乗法等価回路モデルヒストグラムGimp円周率TS-110TS-112PIC16F785LMC662三次元specx.fifortUbuntu疎行列不純物問題Realforceジバニャン方程式ヒストグラム確率論マテリアルデザインP-10境界条件連立一次方程式熱拡散方程式AACircuitHiLAPW両対数グラフ片対数グラフ陰解法MBEナイキスト線図負帰還安定性Crank-Nicolson法EAGLE関数フィッティング

最新コメント
リンク

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