Igorで行列計算

コメントをどうぞ


Igorで行列計算

例として以下の行列T1とT2の掛け算をIgorで計算する。

結果は,

となる。

以上の計算をigorで行った。

以下はそのプログラム例

Make /N=(2,2) T1,T2,T3 //2行2列の行列T1,T2,T3を定義

T1[][0] = {1,2} //0列目に値代入(行列番号は0から始まる)
T1[][1] = {3,4} //1列目に値代入
T2[][0] = {5,6}
T2[][1] = {7,8}

MatrixMultiply T2,T1 //行列の掛け算
T3 = M_product //結果はM_productとして保存される。

AppendToTable T1,T2,T3

matrixmultiplyの計算例

行列計算結果

ただし行列の掛け算部は,以下の演算子”MatrixOp“を用いたほうが簡単。
行列掛け算部の2行は次の1行と等価である。

MatrixOp /O T3 = T2 x T1 //行列の掛け算はx(エックス)を用いる。

MatrixOpを用いる場合には,T3を前に宣言しておく必要なし。
しかし,ここでは逆に,宣言してしまっているから”/O”フラグで上書き許可してる。

Igorで複数軸のグラフ作成

コメントをどうぞ


Igorで複数軸のグラフを作成する方法について記す。

2つの系列でそれぞれ左軸,右軸というグラフはよく作るが,今回は,左軸だけで分離して作成する方法を記す。

 

例として,”Time_s”,”V1″,”V2″という3つのデータ系列を用意し,それをグラフ化する。

“Windows”タブの”New Graph”を選択すると,下図の画面が現れる。”Fewer Choices”ボタンを押すと,簡易設定画面になるが,今回は下図の画面のように詳細設定できる画面ですすめる。

グラフ作成画面

XWaveとしてTime_sを選択し,YwaveとしてV1を選択し,”Add”ボタンを押す。

第2系列を追加

 

同様に,XWaveとしてTime_sを選択し,YwaveとしてV2を選択し,”Add”ボタンを押す。
V2の系列のところで,YAxisの”left”となっているところをクリックし,”New”をクリックする。

そして,Nameの欄に適当に名前をつける。ここでは”L2″とした。

新規軸を追加

 

続いて,”Do It”ボタンを押すと,下図のようなグラフが作成される。

グラフ描画

このままでは,左軸の2つが重なりあってしまい,見栄えが良くないため,修正をする。

左軸をダブルクリック(右クリックして”Axis Properties”でもよい)する。

Axisでleftを選択し,Axisタブの中の”Draw between”のところを0~70%とする。

左軸の修正

 

続いて,Axisで”L2″を選択し,Axisタブの中の”Draw between”のところを80~100%とする。

L2というのは先ほどつけた新規左軸の名前である。

そして,”Distance from Margin”のところを0とする。こうすることで,既存のleft軸と同じ直線上に新規左軸が移動する。

新規左軸の修正

複数軸のグラフはこれで完成である。

あとは,適当にグラフの見栄えを修正すると,下図のようなグラフとなる。

グラフの完成

Igorで複素数計算 ローパスフィルタの特性計算

コメントをどうぞ


Igorで複素数の計算をする方法について書く。

Igorでのコマンドは以下のリンクに詳しい。

http://www.hulinks.co.jp/support/igor/reference/functions/index_complex.html

 

例題として,図のようなローパスフィルタの特性について計算した。

R=100Ω,C=1μFである。

RとCによるローパスフィルタの回路図

RとCによるローパスフィルタの回路図

入力と出力の電圧比は,次式となる。

入出力の電圧比

入出力の電圧比

 

この絶対値が電圧利得の周波数特性になり,極座標表示における,虚部が位相のずれの周波数特性を表す。

Igorでのプログラム例を示す。

用いた関数は,cmplxr2polarである。

cmplx(実部,虚部)

という書式で,複素数を定義する。

wave1=r2polar(wave2)

という書式で,wave2を極座標変換し,wave1を作成する。wave1の実部がwave2の大きさ,wave1の虚部がwave2の位相である。

また,複素数の大きさについては,

wave3=magsqr(wave4)

というようにすればwave4(複素数)の大きさの二乗がwave3(実数)に代入される。

複素数のwave作成時には,Makeコマンドに/cというフラグを付ければよい。

 

プロシージャ画面で,関数(Lowpass)を作成。

ローパスフィルタ特性を計算する関数の定義

ローパスフィルタ特性を計算する関数の定義

実行すると,

関数実行後の特性グラフ化

関数実行後の特性グラフ化

となる。

電圧利得と位相ずれの周波数特性が計算できた。

コマンドを以下に示す。

Function Lowpass()
    Variable R=100
    Variable C=1e-6
    Make /o /N=100 wave_freq = 0.05*x    //0から0.05刻みでwaveの初期値を作成
    wave_freq = 10^(wave_freq)            //logスケールで周波数を定義
    Make /o /N=100 wave_omega = 2*pi*wave_freq
    Make /o /N=100 /c wave_BoardDiag    //入出力の電圧比(複素数)
    Make /o /N=100 mag_BoardDiag        //電圧利得の周波数特性
    Make /o /N=100 phase_BoardDiag        //位相ずれの周波数特性
    wave_BoardDiag = cmplx(1,wave_omega*C*R)    //1+jωC
    wave_BoardDiag = 1/wave_BoardDiag            //1/(1+jωC)
    mag_BoardDiag = real(r2polar(wave_BoardDiag))    //極座標の実部が電圧利得を表す。
    //mag_BoardDiag = sqrt(magsqr(wave_BoardDiag))としても電圧利得が計算できる
    phase_BoardDiag = 180/pi*imag(r2polar(wave_BoardDiag))    //極座標の虚部が位相ずれを表す。
    AppendToTable wave_freq, wave_BoardDiag, mag_BoardDiag,phase_BoardDiag    //各waveの値をテーブルに表示
    Display mag_BoardDiag vs wave_freq                    //グラフ作成
    ModifyGraph grid=1,log=1,tick=2,mirror=1   
    ModifyGraph width={Aspect,1.2}                       
    Label left “V\\Bout\\M/V\\Bin”
    Label bottom “Frequency [Hz]”
    Display phase_BoardDiag vs wave_freq     //グラフ作成
    ModifyGraph grid=1,log(bottom)=1,tick=2,mirror=1
    ModifyGraph width={Aspect,1.2}
    Label left “Phase [deg]”
    Label bottom “Frequency [Hz]”
End

 

Igorで数値積分

コメントをどうぞ


本稿では,Igorで数値積分を実行する。

IgorにはIntegrate1dというコマンドがあり,シンプソン法やルンベルク法を用いた数値積分が簡単に行えるようになっている。
詳しくは,以下のリンクを参照のこと。
http://www.hulinks.co.jp/support/igor/reference/functions/integrate1D.html

例題として,sin(x)の数値積分を行った。
積分範囲は0からπまでとする。
結果は数値積分するまでもないが,2となる。

プロシージャーウィンドウで以下の関数(do_integとinteg_func)を記述する。
関数の名前は任意。
そして,エクスペリメントプロシージャウィンドウで,
Print do_integ(0,pi)
と記述する。
すると答えが得られて,
2
と出力される。

sin(x)の数値積分

 

Function do_integ(xmin,xmax)
    Variable xmin,xmax //積分範囲
    return Integrate1d(integ_func,xmin,xmax,1)  //引数の1はRomberg法の意
End

Function integ_func(inY) //積分したい関数を記述
    Variable inY
 Variable func
 func = sin(inY) //積分したい関数
    return (func) //Integrate1dに値を返す
End

積分範囲を0から2πにしたければ,
Print do_integ(0,2*pi)
とすればよい。
結果はNaN(0の意)が出力される。

積分したい関数をx^2にしたければ,
integ_func関数内のfunc = sin(inY)のところを,
func = inY*inY
とすればよい。

また,以下のように,関数内で積分結果を利用することも当然できる。
先の関数にmain関数を追加して,エクスペリメントプロシージャで,
main()
コマンドを入力すると,
answer = 2
と出力される。

数値積分の他の関数での利用

数値積分の他の関数での利用

Function main()
    Variable imin, imax, c
    imin = 0
    imax = pi
    c = do_integ(imin, imax)
    Printf “answer = %g”, c
End

Function do_integ(xmin,xmax)
    Variable xmin,xmax //積分範囲
    return Integrate1d(integ_func,xmin,xmax,1)  //引数の1はRomberg法の意
End

Function integ_func(inY) //積分したい関数を記述
    Variable inY
 Variable func
 func = sin(inY) //積分したい関数
    return (func) //Integrate1dに値を返す
End

Igor : 拡散方程式 Diffusion equation を解く

コメントをどうぞ


「みその計算物理学」では,一次元の拡散方程式(Diffusion equation)を,
CプログラムやJavaプログラムでの例が紹介されている。
http://www.geocities.jp/supermisosan/kakusan.html

本稿では,Igorで同じ問題を解く(間違っている可能性もあるが)。

一定の温度に熱せられた1次元の系が、両端が温度0にふれているという条件で解く。

Diffusion equation procedure window

拡散方程式(Diffusion equation)のプロシージャ

Function Diffusion_1D_main()
Variable i,j
for (i=0;i<1000;i+=1)
Make /O/N=100 $("wave"+num2str(i))=100 //初期値として100℃にする。
Bound_routine($("wave"+num2str(i))) //境界条件を適用する関数
endfor
for(i=0;i<999;i+=1) //
for(j=1;j<99;j+=1) // 拡散方程式を解く関数Diffuを呼び出す
Diffu($("wave"+num2str(i+1)),$("wave"+num2str(i)),j)
endfor
if (mod(i,100)==0) //iが100でと割り切れる場合に,Tableに追加
AppendToTable $("wave"+num2str(i))
endif
endfor
AppendToTable $("wave"+num2str(i)) //テーブルにwaveを追加する。
Displaywave0,wave100,wave200,wave300,wave400,wave500,wave600,
wave700,wave800,wave900,wave999 //グラフ化する。
End

Function Bound_routine(w) //境界条件を適用する関数
wave w
w[0]=0 //両端の温度を0℃にする。
w[99]=0
end

Function Diffu(w0,w1,jj) //拡散方程式を解く関数
Variable jj
wave w0,w1
w0[jj]=w1[jj]+0.5*(w1[jj+1]+w1[jj-1]-2*w1[jj])
end

図のように,
プロシージャーウィンドウで,
Diffusion_1D_main()
を実行する。

テーブルに計算結果であるwaveが追加され,
グラフ表示される。

Diffusion equation graph and table

グラフとテーブル

Igor:for文を用いてwaveを複数作成する

コメントをどうぞ


C言語でおなじみのfor文を用いてwaveを複数作成する方法を述べる。
後で参照しやすいようにwave名も工夫する。

新規でプロシージャを作成し,for文の構文を確認する。

以下の関数test()を作成する。

Function test() //testという関数作成
    Variable i //変数定義
    for(i=0;i<5;i+=1) //i=0からi<5まで+1ずつ増加
        print i   //iの値を出力
        print “wave”+num2str(i) //数値変数iを文字列に
    endfor  //for文終了
End    //関数終了

print “wave”+num2str(i)について補足。
num2strというコマンドを利用して,数値変数iを文字列に変換する。
そして”wave”という文字列の後につなげることで,
wave1,wave2,…wave9となるようにし,画面に出力する。

エクスペリメントプロシージャでtest()を実行後の結果は次のようになる。

for文の確認

for文の確認

これを踏まえて,waveの作成を行う。

waveの作成はMakeというコマンドを用いる。
Make /N=10 wave2
とすることで,wave2という名前で10行のwaveが作成される。
ここで,先ほど述べたことを応用して,
Make /N=10 $(“wave”+num2str(i))
とする。
$というのは,“文字列をwaveの参照に変換する”コマンドである。
http://www.hulinks.co.jp/support/igor/programming/p_0313.html

作成したwaveをテーブルに追加する。
AppendToTable /O $(“wave”+num2str(i))

次のようにtestという関数を作成する。
Function test()
   Variable i  //変数定義
   for(i=0;i<5;i+=1) //i=0からi<5まで+1ずつ増加
       Make /N=8 /O $(“wave”+num2str(i))  //waveの作成
       AppendToTable /O $(“wave”+num2str(i)) //waveをテーブルに追加
    endfor   //for文終了
End

実行後の結果は次のようになる。

for文を利用したwaveの作成

for文を利用したwaveの作成

Igor:waveの読み込み2

コメントをどうぞ


前回のIgor:waveの読み込みでは,数値データのみの読み込みを行った。
本稿では,ある列にテキストデータがあるデータを読み込む。
例として,以下の図のcsvファイルを読み込む。

読み込むcsvファイルの例
読み込むcsvファイルの例

 

前回のコマンド
LoadWave /J /E=2
では3列目のテキストデータ列が読み込めない。
フラグとして,
/B=”F=0;F=0;F=-2;”
を追加する。
/B“個別の列の特性を指定する”フラグである。
そして,F=<>により,列のデータフォーマットを指定する。

-2 テキスト。
-1 フォーマット不明。Igor がフォーマットを推論。
0 から 5 数値
6 日付
7 時刻
8 日付/時刻
9 8進数
10 16 進数

それぞれの列はセミコロン“;”で区切る。

またFの他に
N=” “とすることで列名を指定できる。
同一列での特性指定の追加する際にはカンマで区切る。
LoadWave /J /E=2 /B=”F=0;F=0;F=-2,N=Day”

 

コマンド実行後のwaveプレビュー
コマンド実行後のwaveプレビュー

列名は”Day”とした。

waveの読み込み完了
waveの読み込み完了

Older Entries Newer Entries