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

2018年7月 3日 (火)

Scratch(プログラム言語)でExponential Sums(指数和)の図示をしてみる。

John Cookさんのブログで知った話題。
毎日Exponential Sumsを描いている。
こういうやつです。
Exponentialsum2
Scratchなら簡単に描ける。
リンクはこちら。
説明はこんな感じで
Turbo mode(Shift+green flag) is recommended.
f(n) : Sums of exp(2*pi*i*(n/dd+n^2/mm+n^3/yy)  , n=1 to nmax
Please try dd=11,mm=21,yy=31 or dd=10, mm=7, yy=1
いろんな絵が描けます。
Exponentialsum_2

2018年6月13日 (水)

Powerpointで数式入力をLaTeX形式で行う。WordはデフォルトでできるけどPowerpointはちょっと初期設定が必要。

Wordはだいぶ前からLaTeX形式で入力できるようになっているのだが、

Powerpointはちょっと初期設定が必要。
制御文字を使わないといけないのだ。
ちょっと説明がわかりにくけど、数式のデザインのツール、の右下をクリックして、
Powerpoint_latex_2
数式オートコレクト にTeX入力の特殊文字を設定する。(左に\TeX, 右に24C9
を入力してからALT+X)
これで、数式入力時に\Tex + スペースでTeXモードになる。
ではやってみよう。
こんな感じで普通にかける(\beginや\endとするとかってに【】になる。\fracとすると〼マークになる)
Powerpoint_latex2

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年9月
            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            
フォト
無料ブログはココログ