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 数値積分 重積分 

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

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAOPアンプPICecalj状態密度常微分方程式モンテカルロ解析トランジスタodeDOSインターフェースPDS5022定電流スイッチング回路分散関係半導体確率論シェルスクリプトレベルシフト乱数HP6632Aバンド構造温度解析可変抵抗I2CR6452AブレッドボードPWscfバンドギャップトランジスタ技術セミナー数値積分反強磁性絶縁熱設計偏微分方程式非線形方程式ソルバフォトカプラシュミットトリガA/DコンバータLM358LED三端子レギュレータカオスISO-I2CGW近似順列・組み合わせマフィンティン半径発振回路74HC4053カレントミラーアナログスイッチPC817CTL431直流動作点解析USBサーボ数値微分固体電子論基本並進ベクトルQuantumESPRESSOスピン軌道相互作用相対論補間FFT開発環境チョッパアンプ単振り子電子負荷BSch量子力学アセンブラパラメトリック解析標準ロジック2ちゃんねるトレーナーバトルポケモンGOスーパーリーグbzqlty状態方程式コバルトLDAイジング模型VESTAブラべ格子六方最密充填構造位相図繰り返し熱伝導Maximaキュリー温度TLP621FETラプラス方程式抵抗失敗談SMPスイッチト・キャパシタgfortranスレーターポーリング曲線ewidth最適化Quantum_ESPRESSOVCAGGA仮想結晶近似ランダムウォークcygwin不規則合金QSGWダイヤモンドQNAPデータロガーマントルシュレディンガー方程式固有値問題自動計測条件分岐格子比熱井戸型ポテンシャルUPSFXA-7020ZRTLP521LM555NE555TLP552Writer509テスタMCU詰め回路三角波過渡解析ガイガー管OpenMPZnOフォノンハーフメタルubuntu最小値最大値熱力学xcrysdenCIF結晶磁気異方性ゼーベック係数スーパーセルUbuntu第一原理計算平均場近似起電力差し込みグラフfsolveブラウン運動フェルミ面awk疎行列縮退文字列不規則局所モーメント入出力化学反応クーロン散乱ヒストグラムRealforceマンデルブロ集合フラクタルキーボード三次元凡例HiLAPW両対数グラフ関数フィッティング熱拡散方程式陰解法片対数グラフグラフの分割シンボルGimp線種Crank-Nicolson法軸ラベル円周率MAS830L負帰還安定性ナイキスト線図EAGLEMBEAACircuitP-10フィルタノコギリ波CapSense2SC1815OPA2277PGALMC662PIC16F785TS-112TS-110等価回路モデルパラメータ・モデル日本語Excel直流解析CK1026トランス連立一次方程式トラックボールPC最小二乗法fcccif2cellPWgui状態図ハイパーリーグPvP擬ポテンシャル不純物問題SIC二相共存重積分電荷密度磁気モーメントEMTO-CPASPR-KKRBroydenTchebyshevルチル構造gnuplot面心立方構造pmixPbTeLmtARTFPLOAMULETsedOpenPBSウルツ鉱構造BaOinterp1初期値ウィグナーザイツ胞L10構造非線型方程式ソルバ正規分布等高線specx.fifortマテリアルデザインヒストグラム確率論ジバニャン方程式TeXFSMedeltquantumESPRESSOリジッドバンド模型スワップ領域岩塩構造デバイ模型半金属全エネルギー固定スピンモーメントc/amultiplot合金境界条件

最新コメント
リンク

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