パソコン・インターネット

2018年5月18日 (金)

円柱導体の高周波インピーダンス(表皮効果含む)の計算をカシオの高精度計算サイトkeisan.casio.jpに自作式として作った。

ちょっと必要になって久しぶりに高周波関係の自作式を作ってみた。
実はベッセル関数の複素引数版が必要で、カシオの高精度計算サイトkeisan.casio.jpにはデフォルトで入っている。便利!
この辺が参考になります。

Skineffect

表皮効果でだんだん抵抗が上がるのと、
内部インダクタンスと外部インダクタンスの和が全インダクタンスになるのでちょっと下がる。
で計算はこちら。

円柱導体の高周波インピーダンス(表皮効果含む)の計算

こんな感じでグラフになる。
94402813825

2018年5月 8日 (火)

Scratch(プログラム言語)でCoupled Map Lattice(CML)を計算してみる。

最近、いろいろGIFアニメをつくるのをよくやっている。
簡単でかつアニメにしておもしろそうなのは、、、結合写像系(Coupled Map Lattice)なんてどうでしょう。
こちらを参考に。
昔、パソコンのスペックが低かった時はこういうのがぴったりのおもちゃだったのでした。
ではスクラッチのリンクはこちら。

Cml

こんな感じでaとgをスライダーで動かせる。ターボモード必須です。
Cmlanime

2018年2月19日 (月)

ソフトバンクが1.7GHz帯を3Gで使えなくしたので、持っているポケットWiFi GL01Pが使えなくなり、代わりに607HWが送られてきた。

ソフトバンクが

一部3Gサービス(1.5GHz帯/1.7GHz帯)提供終了について

としたので、ずっと使っていたPocket WiFi GL01Pが使えなくなった。
どうするかな、と思っていたらYmobileは代替機を送ってきてくれた。

http://www.ymobile.jp/support/relief/nwinfo/1700mhz/

607HWだ(2.1GHz帯を使っていると思う)。
https://wimaxgogo.com/pocket-wifi-607hw/

20180217_121711

これ、最初、何やっても圏外だったのでどういうこと?と思ったが、一枚ペラの移行説明書を見逃していた。その通り設定したら動きました。

2018年2月13日 (火)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その7:メルセンヌツイスタ

高速メルセンヌツイスタについてはこちらを参照:本当に便利に使わせてもらっています。
で、この中のsfmt.zipを使います。この中の
libSFMT.dllをパスが通った場所に置き、
下のSMFTLib.basをインポートして、
    Dim i As Integer
    Dim a As Double
 
    Call InitMt(1)
 
    For i = 1 To 100
        a = NextNormal
        Worksheets("Sheet3").Cells(i + 2, 3) = i
        Worksheets("Sheet3").Cells(i + 2, 4) = a
    Next i
のような感じで使う。
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

2018年2月12日 (月)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その6:ルンゲクッタ8次 DOP853

今回は第6回。ルンゲクッタ8次のDormand & Prince法の有名なFortranのルーチンをVBAに移植しています。
使い方は
DOP853Lib.bas
DriverDOP853Lib.bas
という2つのファイルをインポートします。
DOP853Libのほうが本体で、DriverDOP853Libというのはドライバで問題によって書き換えることになります。
関数はドライバの中のFCNに書きます。
例えばローレンツ方程式なら、
Sub FCN(N As Long, x As Double, y() As Double, F() As Double, RPAR() As Double, IPAR() As Long)
        Dim sigma As Double, r As Double, b As Double
 
        sigma = RPAR(1)
        r = RPAR(2)
        b = RPAR(3)
        F(1) = -sigma * (y(1) - y(2))
        F(2) = -y(2) - y(1) * y(3) + r * y(1)
        F(3) = y(1) * y(2) - b * y(3)
End Sub
密出力ルーチンはSOLOUTで、ここを書き換えて自由なところに出力できます。
Sub SOLOUT(NR As Long, XOLD As Double, x As Double, y() As Double, _
            N As Long, CON() As Double, ICOMP() As Long, ND As Long, _
            RPAR() As Double, IPAR() As Long, IRTRN As Long, XOUT As Double)
' --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS
' --- BY USING "CONTD8", THE CONTINUOUS COLLOCATION SOLUTION
        Dim K As Long
        If (NR = 1) Then
            Worksheets("Sheet2").Cells(JJ + 1, 3) = x
            For K = 1 To N
                Worksheets("Sheet2").Cells(JJ + 1, 3 + K) = y(K)
            Next K
            XOUT = HH
            JJ = JJ + 1
        Else
Label10:
           If (x >= XOUT) Then
              Worksheets("Sheet2").Cells(JJ + 1, 3) = XOUT
              For K = 1 To N
                Worksheets("Sheet2").Cells(JJ + 1, 3 + K) = CONTD8(K, XOUT, CON, ICOMP, ND)
              Next K
              XOUT = CDbl(JJ) * HH
              JJ = JJ + 1
              GoTo Label10
           End If
        End If
End Sub
ドライバ本体の誤差に関わるパラメータは、本家のDOP853の説明を参照。
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

2018年2月 9日 (金)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その5:フィッティング(非線形)

今回の第5回は非線形関数のフィッティングです。
Levenberg-Marquardt Methodを用います。
サブルーチンの引数は前回とよく似ているので同じところは省略ですが、
Sub mrqmin(x() As Double, y() As Double, sig() As Double, ndata As Integer, a() As Double, ia() As Integer, _
      ma As Integer, covar() As Double, alpha() As Double, ByRef chisq As Double, ByRef alamda As Double)
まず初期値を適切に設定し、alamda<0として呼び出します。これで初期化されます。
ここから繰り返しますが、うまく行っていればchi2は小さく、alamdaは10のオーダーで小さくなります。
逆にうまく行ってなければalamdaは10のオーダーで大きくなります。
ここをうまく判定して、最終的に収束したと判定できれば
alamda=0として最後に呼び出し、終了です。
フィッティングするための関数は、VBAが関数を引数に取れないため、名前をきめうちしたサブルーチンとして実装されています。
Sub mrqfuncs(x As Double, a() As Double, ByRef y As Double, dyda() As Double, na As Integer)
となります。パラメータを含む関数形と、そのパラメータごとの微分(ヤコビアン)が必要となります。
例えば、
Sub mrqfuncs(x As Double, a() As Double, ByRef y As Double, dyda() As Double, na As Integer)
 
    Dim i As Integer
    Dim fac As Double, ex As Double, arg As Double
    y = 0#
    For i = 1 To na - 1 Step 3
        arg = (x - a(i + 1)) / a(i + 2)
        ex = Exp(-arg * arg)
        fac = a(i) * ex * 2# * arg
        y = y + a(i) * ex
        dyda(i) = ex
        dyda(i + 1) = fac / a(i + 2)
        dyda(i + 2) = fac * arg / a(i + 2)
    Next i
End Sub
ではガウシアンでフィッティングします。
例題:
(見てもらうと分かるように単に繰り返しているだけで、収束判定はあえていれていません。ここがポイントなので自らのアルゴリズムを組み込んでください)
    Dim n As Integer, m As Integer
    Dim i As Integer, j As Integer
    Dim x(101) As Double
    Dim y(101) As Double
    Dim sig(101) As Double, dyda(7) As Double
    Dim a(7) As Double, ia(7) As Integer
    Dim covar(7, 7) As Double, chisq As Double, ochisq As Double
    Dim alpha(7, 7) As Double, alamda As Double
 
    n = 101
    m = 3
 
 
    For i = 1 To n
        x(i) = Worksheets("Sheet1").Cells(i + 2, 23)
        y(i) = Worksheets("Sheet1").Cells(i + 2, 24)
        sig(i) = 1#
    Next i
 
    For i = 1 To m
        a(i) = 3
        ia(i) = 1
    Next i
 
    alamda = -1
    Call mrqmin(x, y, sig, n, a, ia(), m, covar, alpha, chisq, alamda)
    For i = 1 To 10
        Call mrqmin(x, y, sig, n, a, ia(), m, covar, alpha, chisq, alamda)
    Next i
 
    alamda = 0
    Call mrqmin(x, y, sig, n, a, ia(), m, covar, alpha, chisq, alamda)
 
    Worksheets("Sheet1").Cells(5, 27) = a(1)
    Worksheets("Sheet1").Cells(6, 27) = a(2)
    Worksheets("Sheet1").Cells(7, 27) = a(3)
    For i = 1 To n
 
        Call mrqfuncs(x(i), a, y(i), dyda, m)
        Worksheets("Sheet1").Cells(i + 2, 25) = y(i)
    Next i
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

2018年2月 8日 (木)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その4:フィッティング(線形)

今回の4回目ははフィッティング関数.
(1) 線形 y=a+b*x
(2) 線形    y=a1*f1(x)+a2*f2(x)+・・・+an*fn(x)
(3)  非線形  y=f(a1,a2,・・・, an)
の3種類があります。
まずは線形の(1)はfitというサブルーチンを使います。引数は以下の通り。
Sub fit(x() As Double, y() As Double, ndata As Integer, sig() As Double, mwt As Integer, ByRef a As Double, ByRef b As Double, _
                ByRef siga As Double, ByRef sigb As Double, ByRef chi2 As Double, ByRef q As Double)
データはx(), y()という配列に入っており(全部で ndata個)、
そのデータの標準偏差(yのほう)がsig()に入っているとします。
※ただし、標準偏差がわからないときはmwt=0とすればこのデータを使わずに計算できます。
そのとき、y=a+b*xの形の線形関数にΧ^2を最小化することでフィッティングします。
※パラメータa,bの推定誤差はsiga,sigbに、Χ^2の値はchi2に、当てはまりの良さはqにそれぞれ返ります。
 詳細はNumerical Recipesを参照
事例:
    Dim i As Integer, n As Integer, j As Integer
    Dim x() As Double, y() As Double
    Dim sig() As Double
    Dim a As Double, b As Double
    Dim siga As Double, sigb As Double
    Dim chi2 As Double, q As Double
 
    n = 7
    ReDim x(n), y(n), sig(n)
 
    For i = 1 To 7
        x(i) = Worksheets("Sheet1").Cells(i + 1, 6)
        y(i) = Worksheets("Sheet1").Cells(i + 1, 7)
    Next i
    MsgBox x(1)
 
 
    Call fit(x, y, n, sig, 0, a, b, siga, sigb, chi2, q)
 
    Debug.Print a
    Debug.Print b
 
--
次は一般の線形関数のフィッティングです。
わかりやすいのはこういう多項式フィッティング。
y(x) = a1 + a2*x + a3*x^2 + ・ ・ ・ + aM*x^(M?1)
ですが、より一般的な
y(x)=a1*X1(x)+a2*X2(x)+・・・+aM*XM(x)
の形が可能です。
サブルーチンの引数はこのようになります。
Sub lfit(x() As Double, y() As Double, sig() As Double, ndat As Integer, a() As Double, ia() As Integer, _
    ma As Integer, ByRef covar() As Double, ByRef chisq As Double)
x() ,y()はそれぞれデータが入った配列、Sigはデータの標準偏差で、
パラメータはa()という配列に入ります(ma個)。y=Σa(i)*fi(x)のような形になります。
ただし、ia()という配列を使って0ならそのパラメータはフィッティングしない、1ならする、のように選ぶことができます。
その他パラメータは共分散行列とカイ二乗の値です。Numerical Recipesを参照&あとの事例を参照。
ただし、VBAのサブルーチンが関数を引数に取れない関係上、
フィッティングに使う関数の名前は指定します。ライブラリ中のfuncsというパラメータを
一括で返すサブルーチンを書き換えてください。
多項式の場合は
Sub funcs(x As Double, p() As Double, np As Integer)
    Dim j As Integer
    p(1) = 1#
    For j = 2 To np
        p(j) = p(j - 1) * x
    Next j
End Sub
のようになります。
奇数次のsinで展開するときは、
Sub funcs(x As Double, p() As Double, np As Integer)
    Dim j As Integer
    Dim pi As Double
 
    pi = 3.14159265358979
 
    For j = 1 To np
        p(j) = Sin(x * pi * (2 * CDbl(j) - 1))
    Next j
End Sub
とできる。
矩形波を展開する事例:
    Dim i As Integer, j As Integer
    Dim n  As Integer, m As Integer
    Dim x() As Double, y() As Double
    Dim sig() As Double
    Dim a() As Double, ia() As Integer
    Dim covar() As Double
    Dim chisq As Double
    Dim p() As Double
    n = 21
    m = 4
    ReDim x(n), y(n), sig(n)
    ReDim a(m), ia(m), p(m)
    ReDim covar(m, m)
    For i = 1 To n
        x(i) = Worksheets("Sheet1").Cells(i + 14, 6)
        y(i) = Worksheets("Sheet1").Cells(i + 14, 7)
        sig(i) = 1#
    Next i
 
    For i = 1 To m
        a(i) = 0.1
        ia(i) = 1
    Next i
 
    Call lfit(x, y, sig, n, a, ia, m, covar, chisq)
 
    For j = 1 To m
        Worksheets("Sheet1").Cells(14 + j, 10) = a(j)
    Next j
 
    For i = 1 To n
        Call funcs(x(i), p, m)
        y(i) = 0#
        For j = 1 To m
            y(i) = y(i) + p(j) * a(j)
        Next j
        Worksheets("Sheet1").Cells(i + 14, 8) = y(i)
    Next i
 
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

2018年2月 7日 (水)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数)

今回は3回目でガウス・ジョルダンの消去法(複素数)。
Aをnxn行列、bをnxm行列として、
Ax=bとなるxを計算します。
そして計算したあとのAはもとのAの逆行列で上書きされ、
bは答えのxで上書きされます。
Sub gaussj(ByRef a() As Double, N As Integer, ByRef b() As Double, M As Integer)
Sub Cgaussj(ByRef a() As Complex, N As Integer, ByRef b() As Complex, M As Integer)
実数の場合は
    Dim n as Integer
    Dim Amat() as Double, Bvec() as Double
    n = 4
    ReDim Amat(n, n)
    ReDim Bvec(n, 1)
 
    Amat(1, 1) = 1
    Amat(2, 1) = 1
    Amat(3, 1) = 1
    Amat(4, 1) = 1
 
    Amat(1, 2) = 1
    Amat(2, 2) = 1
    Amat(3, 2) = 1
    Amat(4, 2) = -1
 
    Amat(1, 3) = 1
    Amat(2, 3) = 1
    Amat(3, 3) = -1
    Amat(4, 3) = 1
 
    Amat(1, 4) = 1
    Amat(2, 4) = -1
    Amat(3, 4) = 1
    Amat(4, 4) = 1
 
    Bvec(1, 1) = 0
    Bvec(2, 1) = 4
    Bvec(3, 1) = -4
    Bvec(4, 1) = 2
 
    Call gaussj(Amat, n, Bvec, 1)
 
    For i = 1 To n
        Debug.Print Bvec(i, 1)
    Next i
のように計算できます。
複素数の場合は、
    ReDim A(n, n)
    ReDim b(n, 1)
 
    A(1, 1) = ToComplex(1, 0)
    A(2, 1) = ToComplex(1, 0)
    A(3, 1) = ToComplex(1, 0)
    A(4, 1) = ToComplex(1, 0)
 
    A(1, 2) = ToComplex(1, 0)
    A(2, 2) = ToComplex(1, 0)
    A(3, 2) = ToComplex(1, 0)
    A(4, 2) = ToComplex(-1, 0)
 
    A(1, 3) = ToComplex(1, 0)
    A(2, 3) = ToComplex(1, 0)
    A(3, 3) = ToComplex(-1, 0)
    A(4, 3) = ToComplex(1, 0)
 
    A(1, 4) = ToComplex(1, 0)
    A(2, 4) = ToComplex(-1, 0)
    A(3, 4) = ToComplex(1, 0)
    A(4, 4) = ToComplex(1, 0)
 
    b(1, 1) = ToComplex(0, 0)
    b(2, 1) = ToComplex(4, 2)
    b(3, 1) = ToComplex(-4, -2)
    b(4, 1) = ToComplex(2, 1)
    Call Cgaussj(A, CInt(n), b, 1)
 
    For i = 1 To n
        Debug.Print b(i, 1).x & "+i" & b(i, 1).y
    Next i
 
    For i = 1 To n
        For j = 1 To n
            Debug.Print A(i, j).x & "+i" & A(i, j).y
        Next j
    Next i
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

2018年2月 6日 (火)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その2: FFTライブラリ

昨日の続き。今回はFFT。
Sub FFT(cdata() As Complex, nn As Long, isign As Integer)
のようにサブルーチンとして宣言されています。
まずは複素数のデータ配列を準備します。2のべき乗個(Longで宣言)のデータが必要です。
Dim z1() as Complex
Dim n as Long
n=4    ' 2のべき乗
Redim z1(n)
のように宣言して、データをいれ、
Call FFT(z1, n, 1)
のように呼び出します。z1にデータが上書きされます。
最後の1を-1にすると逆変換になりますが、
データ個数nでは割られていません!必要ならば後で自ら割る必要があります。
(これはNumerical Recipesのルーチンの仕様です)
なので変換→逆変換を繰り返すと元に戻らずにn倍された値になります。
例題:
    z1(1) = ToComplex(1#, 2#)
    z1(2) = ToComplex(2#, 3#)
    z1(3) = ToComplex(3#, 4#)
    z1(4) = ToComplex(4#, 5#)
 
    Call FFT(z1, n, 1)
    Call FFT(z1, n, -1)
 
    For i = 1 To n
        z1(i) = Cdiv(z1(i), ToComplex(CDbl(n), 0#))
    Next i
    For i = 1 To n
        Debug.Print z1(i).x & "+i" & z1(i).y
    Next i
とすると元に戻ります。
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

2018年2月 5日 (月)

Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その1:複素数演算の基礎

会社では、まあ開発部門は何のプログラム言語を使ってもいいのだけれど
工場ではそういうわけにもいかない。そもそも外部のネットワークに接続できなかったりする。
ちょっとした設計や計測に使うソフトウェアを作って配りたい、
というときに非常に困る。
ソフトウェアを買え、という話もあるが、損益が厳しい工場が買うはずがないのだ。(これを導入したらいくらの売り上げUPか利益率いくらUPが定量的に示せないと何も買わないのだ)
MicrosoftのOfficeくらいならどの工場のPCでも入っている、ということで
フルスクラッチでExcelのマクロ(VBA)で計算・計測プログラムを作って配っていたりする。
DLLですらダウンロードできないので、全部ソースに書いているという。
※ちなみに仕事でやっているわけでなくて、単なるサービスというか、趣味。
 私はもはやエンジニアでもないのだったりする。
ということで、今回はそのライブラリを紹介します。まあそれなりにみんなに使ってもらっているので致命的なバグはないと思いますが、細かいミスはご勘弁を、、、(責任は全くとらない方向で)
まずライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用
――
今回は第一回ということで複素数演算などの使い方。
Numerical Recipesについてはこちらから本が丸ごと読めます。
まずはVBAの編集画面から、ファイルのExportを行い、
ComplexMath.bas
を読み込みます。これで複素数演算(FFT,連立一次方程式,逆行列)、フィッティング(線形・非線形)などが
使える状態になります。
複素数は構造体で定義されています。
Public Type Complex
    x As Double
    y As Double
End Type
この変数を作るときは、
Dim x As Complex
のように宣言します。
もちろん、配列も宣言できます。行列を作るには
Dim A() as Complex
Redim A(3,3)
とすれば可変長の行列ができます。
※Numerical Recipesの表現に合わせたので、0からではなくて1からです。
 明示的に最初に Option BASE 1 としておくと便利
値を代入するには、
z = ToComplex(1.2, 3.5)
のように実部と虚部に分けて設定することができます。
この場合 z=1.2+3.5i と入力したことに相当します。
また
z = x
のように別の複素数は直接代入できます。
四則演算については +,-,×,÷はそれぞれ
    z = Cadd(x, y)
    z = Csub(x, y)
    z = Cmul(x, y)
    z = Cdiv(x, y)
Cdiv(x,y)のみ順番に注意が必要で、x/yの意味です。
複素共役,絶対値、偏角はそれぞれ、
z=Conj(z), a=Cabs(z), p=CPhase(z)
実数部分、虚数部分を取り出すにはそれぞれ
a=Rez(z), b=Imz(z)
などとします。
※ついでに角度が±πに対応したatan2(x,y)という関数も使えます。x/y=tanθの逆関数。
複素数の指数、対数関数、べき乗関数はそれぞれ(多価関数は±πのみを取る)
z=Cexp(x), z=Clog(x), z=Cpow(x,y)となる。ただしべきはz=x^yの意味。
--
では次回はFFTを。
--
こちらも参照:

より以前の記事一覧

2018年5月
    1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31    
フォト
無料ブログはココログ