Scilabで繰り返し計算: ロジスティック写像

Scilabを使えば、非常に簡単に数値計算のプログラムを書くことができます。例えば、たった8行のプログラムでFig.1のような計算が出来てしまいます。

今回は、プログラミングの基本のひとつであるループ処理の例としてロジスティック写像の計算を行いました。Scilabはインタープリタなので、ループ処理は苦手です。そのため出来る限りループを無くすようにする実例にもなっています。

001_20130520003209.png
Fig.1: ロジスティック写像。横軸が繁殖率a、縦軸が100世代目から200世代目までの個体数X(n)。



ロジスティック写像


生物の個体数Xが世代nに従ってどのように変化するかという研究からロジスティック方程式というものが知られています。

X(n+1) = a \times X(n) \times (1-X(n))

この式の見た目はとても簡単であるにもかかわらず、個体数の初期値X(1)や繁殖率aに対して、充分時間がたったときの個体数(例えば100世代目の個体数X(100)など)が非常に敏感に反応するため、カオス理論と関連して議論される面白い方程式となっています。(参考:ロジスティック写像:長崎県立大学伊藤研究室)

002_20130520024958.png
Fig.2: 繁殖率aに関して 3 < a < 4 の部分を拡大したもの


また、数値計算の例題としてもシンプルな繰り返しのコードで思いもよらないほど複雑な、かっこいい結果が得られるので良く見かけます。Fig.1,2はロジスティック方程式を繁殖率aを0<a<4の間で計算したもので、縦軸に個体数X(n)を100 < n < 200の範囲でプロットしてあります。個体数の初期値はX(1)=0.2としました。
本エントリではScilabを用いたループ計算の例としてこのロジスティック写像のプログラムを紹介します。特にScilabのようなインタプリタではループ処理と非ループ処理で計算速度に大きな差が出るため、この点についても書きます。

繰り返し計算


まずは、簡単のために繁殖率をa=3.5に固定して、世代数nを順次変化させるプログラムを書きます。
Scilabで繰り返し処理をするにはforを使います。(参考:繰り返し:Scilab簡易リファレンス-PukiWiki)

nmax = 50;

a = 3.5;
x = 0.2 * ones(nmax,1);

for n = 1: nmax - 1 do
x(n+1) = a * x(n) * (1.0 - x(n));
end

plot(x,'o-r');
xlabel("n");
ylabel("X(n)");

plot([0,0],[0,1]);


最後の行の(0,0)と(0,1)のプロットは描画領域を調整するためのダミーです。


003_20130521000718.png
Fig.3: a=3.5 X(1)=0.2のときの個体数X(n)の世代発展


Fig.3を見ると30世代ぐらいまでは過渡領域で、それ以降は4つの値を行き来していることが分かります。
Fig.1-2は充分時間がったったときの個体数で、a=3.5において4つの値を行き来する事は四本の線に分かれている事から読み取れます。

繰り返し計算の入れ子


それでは次にFig.1-2を描くプログラムについて考えます。
先ほどは世代数nを変更するために繰り返し命令forを使いました。このことを応用してforを入れ子にすれば、さらに繁殖率aを変更した繰り返し計算のプログラムがかけます。

nmax = 200;
nout = 100;

A = [0:0.001:4];
X = 0.2 * ones(nmax,length(A));

for m = 1:length(A) do
for n = 1: nmax - 1 do
X(n+1,m) = A(m) * X(n,m) * (1.0 - X(n,m));
end
end

plot(A,X(nout:nmax,:),'.r','markersize',1);

このプログラムは悪くないのですが、ベストではないです。

繰り返し計算は出来れば避ける


ループ計算は遅いので,できるだけ使わないでは、同じ内容の計算であってもForループの有無により数10倍計算速度が変わる例が示されています。

ループを使わずに済ませることができるかどうかはケース・バイ・ケースなのですが、今回は入れ子にせずとも1段階のループで処理することが出来ます。
今回のロジスティック写像のプログラムは、私の環境では数10倍どころか、ループ部分の計算時間が約14秒から約0.03秒まで短縮できました。(もっとも、そのあとのグラフの描画に数秒程度かかってしまうのですが。)

nmax = 200;
nout = 100;

A = [0:0.001:4];
X = 0.2 * ones(nmax,length(A));

for n = 1: nmax - 1 do
X(n+1,:) = A .* X(n,:) .* (1.0 - X(n,:));
end

plot(A,X(nout:nmax,:),'.r','markersize',1);


ただし、Scilabは細かいことを気にせずにプログラムが書けることが最大の利点なので、計算速度向上のために腐心して時間を無駄遣いするぐらいなら、多少効率が落ちようともさっさとコードを書いてしまうほうがいいという面もあります。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 繰り返し カオス 

Scilab入門書籍2冊

CQ出版のツール活用シリーズからScilabの入門書が発売されました。波形解析のための数値計算ソフト Scilab入門: スペクトラム/ノイズ分析から特徴抽出まで (ツール活用シリーズ)です。

内容は実験データの解析で移動平均やフーリエ解析、ウェーブレット解析を扱っています。Scilabもそのうちブログで扱おうと思っており、丁度良い機会のなのでカットシステムから発売されているScilab入門―電気電子工学で学ぶ数値計算ツールもあわせて2冊紹介します。




数値計算ソフト


本格的なプログラミング言語(例えば、C言語やFortran)を使わずに、簡単な数値計算をする場合、どのようにするのが良いでしょうか?

表計算ソフトのExcelを使うのもひとつの手です。例えばExcelで操る! ここまでできる科学技術計算では科学技術計算の90%はExcelで対応できる!としています。

実際これはその通りで、大部分の計算はExcelで充分計算可能だと思います。ただ、残りの10%の計算が必要な場合もありますし、なによりExcelで対応可能であっても他のソフトで計算するほうが結果的に簡単だということもしばしばです。

そんなとき、より数値計算に向いたソフトが使われます。

Scilabは無料の数値計算ソフトのひとつです。
有償の数値計算ソフトとしてはMATLABが有名で、Scilabは互換ではないものの、大体似たような事が出来ます。

具体的な応用例として、電子工作関連ならのりたんさんがscilabで遊ぼう (1)の一連のエントリにてPID制御のシミュレーションを行っています。

他にもなどさまざまなことが出来ます。


001_20130519233508.png
Fig.1: Scilabによる三次元プロットのデモ


Scilab関連書籍


Scilabの関連書籍は、例えばAmazonでScilabを検索してみると相当数が見つかりますが、その実は、例えば本の主眼を数値計算においていて、Scilabの使い方については簡単に触れられているだけといったものが多いです。

そんな中で、私がねがてぃぶろぐの電子工作ブログという側面から1冊だけ入門書を挙げるならScilab入門―電気電子工学で学ぶ数値計算ツールです。
この本は、Scilabの基本的な使い方から、様々な数値計算までを実例とソースコードを交えて解説しています。
回路設計前のシミュレーションと同様、モデルから出発して値を計算することに主眼が置かれています。

一方で、今回新しく発売された波形解析のための数値計算ソフト Scilab入門: スペクトラム/ノイズ分析から特徴抽出まで (ツール活用シリーズ)は、既に実験などから得られているデータを解析することに主眼を置いた解説書です。



言うなればScilab入門―電気電子工学で学ぶ数値計算ツールは計算屋向けで波形解析のための数値計算ソフト Scilab入門: スペクトラム/ノイズ分析から特徴抽出まで (ツール活用シリーズ)は実験屋向けでしょうか。いや、違うかも。実験屋さんでも簡単に数値計算できるから前者がオススメなのだし。

参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 

擬似Cカーブ可変抵抗の定数設計

LTspiceで擬似Cカーブ可変抵抗コメント欄にて、擬似Cカーブ可変抵抗の定数設計の実際について質問をいただきました。

擬似Cカーブ可変抵抗(Rc)は、Bカーブ可変抵抗(Rb)と並列に入れる固定抵抗(Rp)から構成されるためパラメータが二つ存在します。この二つのパラメータを上手に選ぶことによって擬似Cカーブ可変抵抗の抵抗値(Rc)とその曲率の二つを変化させることが出来ます。

001_20130515013645.png 002_20130515013644.png


擬似Cカーブ可変抵抗


電子回路の「ボリューム・つまみ」として可変抵抗は頻繁に使用されますが、そのつまみの回転角と抵抗値の間の関係は『直線的なBカーブ』『指数関数的なAカーブ』に加えて『対数関数的なCカーブ』の3種類が存在します。

001_20090216001312.png


しかしながら、前者2つと比べてCカーブの可変抵抗はほとんど使われることが無いため、電子工作部品としても入手がやや困難です。

このためLTspiceで擬似Cカーブ可変抵抗では、Bカーブの可変抵抗と並列に普通の固定抵抗を接続することで、抵抗値の変化を上に凸なCカーブ的なものに出来る事をシミュレーションから確認しました。

この記事に対して、コメント欄にて、通りすがりさんに設計に関する質問をいただきました。

はじめまして。
50kCカーブ2連ポットを探していてたどり着きました。
記事を拝見しましたが難しくて(^^;
50kCカーブを作ろうとしたら100kBカーブにどのくらいの抵抗を抱かせればいいのでしょうか?
ご教授いただければ幸いです。


そこで、今回のエントリでは、具体的に擬似Cカーブの可変抵抗を作るに当たっての定数設計に関して書きます。

並列抵抗Rpの計算式


設計する擬似Cカーブ可変抵抗の値をRc、使用するBカーブ可変抵抗の値をRb、並列に入れる固定抵抗の値をRpとします。

擬似Cカーブ可変抵抗Rcは、RbとRpの並列接続なのでその値は
\frac{1}{R_c}=\frac{1}{R_b}+\frac{1}{R_p} ・・・(1)
から
R_c = \frac{R_b R_p}{R_b + R_p} ・・・(2)
であると分かります。

これをRpについて解くと
R_p = \frac{R_b R_c}{R_b - R_c} ・・・(3)
となります。

Bカーブ可変抵抗Rbと並列抵抗Rpの決定


今回の計算例として通りすがりさんの値Rc=50kΩを使います。

(2)式からRcを決めるための未知数はRbとRpの2つあります。しかしながら、Rbは市販のBカーブ可変抵抗なので、おのずと選ぶことの出来る値が限られてきます。
入手性が良さそうで、並列にして50kΩが作れそうな抵抗値として、差し当たりRb = 100kΩ、150kΩ、200kΩ、500kΩあたりを検討することにします。

ここまででRc = 50kΩとRbの候補が4種類決まるので(3)式から、それぞれのRbに対応したRpが計算できます。

Bカーブ抵抗(Rb)並列抵抗(Rp)
100kΩ100kΩ
150kΩ75kΩ
200kΩ66.7kΩ
500kΩ55.6kΩ
table.1: Rc = 50kΩのときのRbとRpの組み合わせ


擬似Cカーブ可変抵抗の曲率


table.1に50kΩの擬似Cカーブ可変抵抗を構成することが出来るBカーブ可変抵抗と並列抵抗の組み合わせの例を挙げました。

では、これらは全て同じ特性となるのでしょうか?
―――答えはNOです。

このことを確認するためにLTspiceで擬似Cカーブ可変抵抗とほぼ同じ方法でLTspiceを用いて擬似Cカーブ可変抵抗の曲率を計算してみました。

ただし、グラフの横軸をBカーブ可変抵抗器の抵抗値にしてしまうとRbの値によって横軸がそろわなくなるので、多少トリッキーではありますがLTspiceで電圧制御抵抗(VCR)の方法を使いました。
fig.2の横軸の単位が電圧となっていますが、これは電圧ではなくBカーブ可変抵抗の回転角であると考えてください。0から1で端から端まで回しきるイメージです。

001_20130515013645.png
fig.1: 擬似Cカーブ可変抵抗のスケマティック

002_20130515013644.png

fig.2: 抵抗の組み合わせによる曲率の違い。横軸は電圧ではなく、可変抵抗の回転角と読み替えてください。(0で抵抗値が最小、1で最大。)


シミュレーション結果から、どの組み合わせであっても、つまみを最小から最大に回しきった時に設計どおり抵抗値が0Ωから50kΩまで変化することがわかります。
しかしながら、その曲率は抵抗の組み合わせによって異なっており、Bカーブ可変抵抗に大きい抵抗値を選ぶほど曲率が大きくなっています。

本物のCカーブ可変抵抗がこの中でどの曲率に一番近いのかは、残念ながら私は知りません。

対数的な変化をするCカーブ可変抵抗とは逆に、指数関数的な変化をするAカーブ可変抵抗は、オーディオの音量変化に頻繁に利用されます。これは、人間の五感が対数的な挙動を示すことが理由です。つまりオーディオのボリュームつまみに指数関数的なAカーブ可変抵抗を利用すると実際の音量も指数関数的に増えているにもかかわらず、人間の聴感からは直線的に音量が増えたように感じるということです。

Cカーブ可変抵抗もまた同様に、人間がつまみを回したときに『直線的に変化したな』と感じるような曲率になっているのが望ましいはずです。
実際のところ、私は擬似Cカーブ抵抗の使い道を自分でも良くわかっていないのですが、どういった曲率が人間にとって一番『しっくりくる』のかはカットアンドエラーで決めても良いのではないかと思います。

関連エントリ




付録


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


フィードバック



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

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


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


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

tag: LTspice 可変抵抗 

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関数フィッティング

最新コメント
リンク

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