Scilabの論理演算で条件分岐

Scilabで条件分岐:テント写像ではif then elseを利用して条件分岐のプログラムを書きました。
しかしながらScilabで繰り返し計算: ロジスティック写像で指摘したとおりScilabではループ計算をさせるととたんに計算速度が遅くなります。

前回のプログラムではif then elseを使うためにforループを二重になってしまっています。
そこで今回は、論理演算を用いてループが減るようにしました。

texclip20130718112739.png

nmax = 200;
nout = 100;

Mu = [1:0.001:2];
X = 0.5 * ones(nmax,length(Mu));

for n = 1: nmax - 1 do
X(n+1,:) = Mu .* X(n,:) .* (X(n,:) < 0.5) + ..
Mu .* (1.0 - X(n,:)) .* (X(n,:) >= 0.5);
end

plot(Mu,X(nout:nmax,:),'.r','markersize',1)
xlabel("$\mu$");
ylabel("$x_n$");



論理演算


Scilabでは結果をT(true, 真)とF(false, 偽)で返す論理演算があります。
例えば、以下のようなものです。
M = [0:0.1:1]
M > 0.5
この結果は F F F F F F T T T T T といった様に要素にTとFを持つ行列になります。

これを F=0,T=1 の数値としてそのまま計算に使うことが出来ます。例えば
M = [0:0.1:1]
1 * (M > 0.5)
の結果は 0 0 0 0 0 0 1 1 1 1 1 です。

関連エントリ




参考URL




付録


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




参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 条件分岐 繰り返し カオス 

Gimpで回転する地球

GIMPで回転する地球のgifアニメーションを作成する方法・前編後編の通りに回転するgifアニメーションを作ってみました。カッコイイ!

globe.gif

Fig.1 自転する地球


関連エントリ




参考URL




フィードバック



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

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


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


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

tag: Gimp 

おまえは今まで食ったパンの枚数をおぼえているのか?

購入する食パンは常に6枚切りとし、バターは8個入りである。食パン1枚につき常に1個のバターを塗ることとし、バターは食パンに塗る以外の用途で消費することはないとする。

食べた数食パンの残りバターの残り
068
157
246
335
424
513
662
751
848
937
1026
1115
1264
1353
1442
1531
1628
1717
1866
1955
2044
2133
2222
2311
2468
table.1: パンの枚数とパンとバターの残量


と言ったようにすれば、パンの枚数を覚えていなくてもある程度の制約が可能かと思ったのですが、実際に冷蔵庫を見てみるとなぜか偶数枚のパンと奇数個のバターが・・・一体何があったのか?どこかの棚の下に溶けたバターが「べちゃー」っとなってるとかだと嫌だなぁ。

参考文献


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

たしかなまんぞく

間々あることなのですが急に漫画が読みたくなったので、ワールドトリガー(1-5巻 以降続巻)を購入しました。





なかなか面白いです。
ジャンプ王道のバトル漫画ですが、主人公(メガネ君の方)は、はっきり言って弱っちいので良くあるインフレバトルにはならなそうなのも好印象。
まだ全然ストーリが始まったばかりなので、続きが楽しみです。

Scilabで数値微分 その2

Scilabで数値微分 その1に引き続き数値計算の常識に従ってf(x)=sin(x)の数値微分を行いました。

今回は刻み幅hを可変させたときのf(x,h)の値の変化から、誤差が最小となるhの最適値を探しました。


Scilabで数値微分 その1では数値計算の常識に従ってf(x)=sin(x)のx=0.3πにおける微分値を計算しf'(x)=cos(x)との比較を行いました。

その結果、微分を差分に置き換える際の刻み幅hには最適値があり、大きすぎても小さすぎても良くないということがわかりました。

Scilabで数値微分 その1はf(x)=sin(x)のような簡単な関数だったので、あらかじめ解析的に微分を行うことが出来ました。しかしながらScilabで金属の電子比熱の様にf(x)の計算自体に非線形方程式ソルバを使うような複雑な場合は、真の値と比較して最適なhを探すことは出来ません。むしろ、解析的に微分が出来るような場合は、そもそも数値的に微分をする必要が無いとも言えるので、解析的に微分できないときにどのように最適なhを探すのかというのは重要な問題です。

実践的な刻み幅の決定法


数値計算の常識には以下の様にあります。

実践的には,f1(x,h)とf1(x,2h)との"差"を観察して,
「hを半分にしたとき"差"が約半分になる傾向が保たれる範囲でもっとも小さなh」
をえらぶのがよい.


そこでこの計算をScilabで行います。

Scilabプログラム


作成したScilabのプログラムはdiff2_sce.txtとなりました。

clear;

// *** 刻み幅の設定 ***
n = [0:1:50];
h = 2 ^ (- n);

x = 0.3 * %pi;

// *** 微分の近似値 ***
// 前進差分
f1 = (sin(x + h) - sin(x)) ./ h;
f12 = (sin(x + 2 * h) - sin(x)) ./ (2 * h);
// 中心差分
f2 = (sin(x + h) - sin(x - h)) ./ (2 * h);
f22 = (sin(x + 2 * h) - sin(x - 2 * h)) ./ (4 * h);
// 前進差分に対するRomberg 1段公式
romberg1 = 2 * ((sin(x + h) - sin(x)) ./ h - 0.5 * (sin(x + 2 * h) - sin(x)) ./ (2 * h));
romberg12 = 2 * ((sin(x + 2 * h) - sin(x)) ./ (2 * h) - 0.5 * (sin(x + 4 * h) - sin(x)) ./ (4 * h));
// 中心差分に対するRomberg 1段
romberg2 = ((sin(x + h) - sin(x - h)) ./ (2 * h) - 0.25 * (sin(x + 2 * h) - sin(x - 2 * h)) ./ (4 * h)) / 0.75;
romberg22 = ((sin(x + 2 * h) - sin(x - 2 * h)) ./ (4 * h) - 0.25 * (sin(x + 4 * h) - sin(x - 4 * h)) ./ (8 * h)) / 0.75;

// *** グラフの軸設定 ***
a = gca();
a.data_bounds = [0,1E-14;50,1];
a.log_flags = "nl";

// *** グラフのプロット ***
plot(n,abs(f1 - f12),'-sr');
plot(n,abs(f2 - f22),'-sm');
plot(n,abs(romberg1 - romberg12),'-sb');
plot(n,abs(romberg2 - romberg22),'-sg');

// *** グラフの体裁 ***
xlabel("n");
ylabel("Err");
legend(['f1(x)','f2(x)','Romberg f1(x)','Romberg f2(x)'],4);


結果


Scilabで数値微分 その1で計算した数値微分と解析解の比較「f(x,h)-f'(x)」をFig.1に、今回計算した「f(x,h)-f(x,2h)」をFig.2に示します。これらは非常に良く似た形をしていますが全くの別物です。

001_20131110042107c73.png

Fig.1: 数値微分と解析解の比較

002_20131110080640da7.png
Fig.2: 数値微分の刻み幅を2倍にしたときの変化の度合い


数値計算の常識には以下の様にあります。

このような"差分"が「規則的に変化している」ことが,「打切り誤差が丸め誤差より優越していて,通常の打切り誤差の漸近理論が適用できる」ことの実際的な判別法として役立つ.


Fig.2をみるとnを大きくしていっても(つまりhを小さくしていっても)規則的に"差分"が小さくなっていっている部分が、丸め誤差の影響を受けず、打切り誤差が支配的となっている領域です。この範囲でもっとも大きなn(もっとも小さなh)を選んだとき、実際の誤差がほぼ最小になっていることがFig.1との比較から読み取れます。

関連エントリ




付録


このエントリで使用した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関数フィッティング

最新コメント
リンク

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