雑感等

音楽,数学,語学,その他に関するメモを記す.

波動方程式の有限差分法を半自動で導出(sympy)

sympyで差分化した2次元の波動方程式を定義して、式の整理を自動で実行する。

jupyterで実行時、最終的な実行結果は下記。

uの添え字は、1番目が時間、2番目と3番目がxとy。

左辺がu_{1,0,0}なので、次の時間ステップにおける波の振幅を表す。

右辺は、u{0,?,?}とu{-1,?,?}で構成されているため、既知の値から計算できる。

jupyter labから.mdで出力したものをそのままここに貼っても、texの数式がうまく表示されなかった。 下ページに従い、pythonコード形式に変換した。 ipynb ファイルを py ファイルに変換する方法 | ふじやまエッグのブリコルール日誌

#!/usr/bin/env python
# coding: utf-8

# # solve
# $$
# \frac{\partial^2 u}{\partial t^2} = c^2 \left[ \frac{\partial^2 u}{\partial x^2}+ \frac{\partial^2 u}{\partial y^2} \right]
# $$
# 
# with explicit finite difference method

# In[1]:


import sympy as sym


# ## define symbols

# discretized wave: $u_{time, x, y}$
# 
# mesh size of space: $ dx $
# 
# mesh size of time: $ dt $
# 
# speed of wave: $c$

# In[2]:


u =[ [[0 for i in range(3)] for j in range(3)]for k in range(3)]

for nt in range(3):
    for nx in range(3):
        for ny in range(3):
#             u[nt][nx][ny] = sym.Symbol("u[nt+("+str(nt-1)+"),nx+("+str(nx-1)+"),ny+("+str(ny-1)+")]") # for python code
            u[nt][nx][ny] = sym.Symbol("u_{"+str(nt-1)+","+str(nx-1)+","+str(ny-1)+"}") # for readable math expression
    pass
u[0][0][0]


# In[3]:


dx = sym.Symbol("dx")  # mesh size of space
dt = sym.Symbol("dt")  # mesh size of time
c = sym.Symbol("c")  # wave speed


# ## define partial derivative
# 
# $$ \frac{\partial^2 u}{\partial x^2} := \frac{-u_{t, x-1, y} + 2u_{t,x,y}-u_{t,x+1,y}}{dx^2} $$
# $$ \frac{\partial^2 u}{\partial y^2} := \frac{-u_{t, x, y-1} + 2u_{t,x,y}-u_{t,x,y+1}}{dx^2} $$
# $$ \frac{\partial^2 u}{\partial t^2} := \frac{-u_{t-1, x, y} + 2u_{t,x,y}-u_{t+1,x,y}}{dt^2} $$
# 

# In[4]:


d2x = sym.expand(((u[1][1][1] - u[1][0][1]) / dx - (u[1][2][1] - u[1][1][1]) / dx) / dx)
d2x


# In[5]:


d2y=sym.expand(((u[1][1][1] - u[1][1][0])/dx-(u[1][1][2] - u[1][1][1])/dx)/dx)
d2y


# In[6]:


d2t=sym.expand(((u[1][1][1] - u[0][1][1])/dt-(u[2][1][1] - u[1][1][1])/dt)/dt)
d2t


# ## define discretized wave equation

# In[7]:


eq=sym.Eq(d2t,c**2*(d2x+d2y))
eq


# ## solve the above eq. for $u_{1,0,0}$=u[2][1][1]

# In[8]:


rhs=sym.solve(eq,u[2][1][1])[0].expand().simplify()
rhs_simp=sym.collect(rhs,c).collect(dt).collect(dx)
rhs_simp


# # discretized wave eq. for computation
# 
# lhs is unknown value
# 
# rhs is constructed with known value

# In[9]:


sym.Eq(u[2][1][1],rhs_simp)

lilypondコンパイルエラーの対処 error: syntax error, unexpected SCM_TOKEN, expecting '='

コンパイルエラーの修正

lilypondファイルをコンパイルすると、以下のエラーメッセージが出力。

././01-byrd-a3-0-score.ly:36:62: error: syntax error, unexpected SCM_TOKEN, expecting '='
                \override StaffGrouper #'staff-staff-spacing
                                                             #'padding = #4.5

"01-byrd-a3-0-score.ly"中の以下の文を書き換えた。

\override StaffGrouper #'staff-staff-spacing #'padding = #4.5

\override StaffGrouper.staff-staff-spacing.padding = #4.5


詳細

コンパイルしたいlilypondファイル:

Mass for 3 Voices (Byrd, William) - IMSLP: Free Sheet Music PDF Download

このページにあるEngraving files (Lilypond)

編集者    Allen Garvin
出版社情報 Dallas: Hawthorne Early Music, 2016.
著作権   
Creative Commons Attribution-NonCommercial 4.0

コンパイルエラー:

ダウンロードしたzipを展開し.\single-parts\01-byrd-a3-0-score.lyコンパイルしたが、下記のエラーメッセージ。

GNU LilyPond 2.20.0
Processing `./01-byrd-a3-0-score.ly'
Parsing...
././01-byrd-a3-0-score.ly:36:62: error: syntax error, unexpected SCM_TOKEN, expecting '='
                \override StaffGrouper #'staff-staff-spacing
                                                             #'padding = #4.5
././01-byrd-a3-0-score.ly:36:62: warning: Ignoring non-music expression
                \override StaffGrouper #'staff-staff-spacing
                                                             #'padding = #4.5
././01-byrd-a3-0-score.ly:36:74: warning: Ignoring non-music expression
                \override StaffGrouper #'staff-staff-spacing #'padding =
                                                                         #4.5

・修正方法の検討:

lilypond公式ドキュメント LilyPond 記譜法リファレンス: 5.3.7 連想配列を変更する のコード例では、

\override StaffGrouper.staff-staff-spacing.basic-distance = #7

のように書いてあり、StaffGrouperstaff-staff-spacing連想配列のキーbasic-distanceの間がピリオドで繋がれている。

・対処:

そのため、\override StaffGrouper #'staff-staff-spacing #'paddingテキストエディタ\override StaffGrouper.staff-staff-spacing.paddingに置換。

置換したファイルをコンパイルすると、コンパイルエラーは出ず、imslpにあるpdfと同様のファイルが生成される。

楽譜を見ながらMass for 3 Voicesを聞こうとしたら調があっておらず、楽譜を移調しようとした。

imslpにlilypondのファイルがあり、ソースを修正して手元でコンパイルしようとしたところ上記エラーが発生。

波動方程式を有限差分法の陰解法で解く(Julia)

波動方程式を↓の文献にある式で数値計算する。

https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ws2016-2017/num_methods_i/wave.pdf

https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/part2_hypebolic/node20.html

記号

  • uは波の変位
  • uの右上の添え字(j-1, j, j+1)は離散時間の番号
  • uの右下の添え字(i-1, i, i+1)は離散空間の番号
  • Δtは離散時間の単位時間(微小時間)
  • Δxは離散空間の単位長さ(微小長さ)
  • cは波の伝播速度

式の意味

文献の式の右辺を少し変形すると↓になる。

f:id:kazmus:20210821145200p:plain

左辺は \frac{\partial^{2} u(x,t)}{\partial t^{2}}を現在時刻(j)で中心差分に置き換えたもの。

右辺は \frac{\partial^{2} u(x,t)}{\partial x^{2}}を、1微小時刻前(j-1)における中心差分に置き換えたものと、1微小時刻後(j+1)における中心差分に置き換えたものの平均。

現在の値を、過去と未来の平均に置き換えられる理由は理解していない。もしかしたら冒頭の文献にあるかもしれない。

(陽解法の場合は右辺は現在時刻(j)だけで表し、 l.h.s. = c^{2} \cdot\frac{u^{j} _{i+1} - 2u^{j} _{i} + u^{j} _{i-1}}{\Delta x^{2}}

行列形式にしやすいように変形

さらに変形して↓になる。

f:id:kazmus:20210821150204p:plain

ここで、定数部分を↓としてまとめて表示する。

f:id:kazmus:20210821150237p:plain

さらに、 u^{j+1}_{??}の項(未来の項・未知数)を左辺に、それ以外を右辺に集めると↓。

f:id:kazmus:20210821152657p:plain

今回の境界条件

 u^{j}_{i}の添え字の範囲は

  • 時刻  j=0 \cdots J
  • 空間  i= 0 \cdots I

境界条件

  • 初速度0: u^{-1}_{i} = u^{0}_{i}
  • 両端固定: u^{j}_{-1} = u^{j}_{I+1}=0

行列表示

境界条件を適用して、行列を使って整理すると

f:id:kazmus:20210821225341p:plain これを使ってJuliaで実装する。

Juliaプログラムソースと実行結果

実行結果↓ 横軸が空間方向、縦軸が時間方向、色が変位の大きさ。

f:id:kazmus:20210821230417p:plain 左右端で位相が反転し、固定端反射している。

図の下方(初期状態付近)だと孤立波の形を保っているが、時間が経過すると波形が崩れる。

ソース↓

using Plots
using WAV
using LinearAlgebra

const t_max = 3000  #時間点数
const x_max = 200   #空間点数

const dt = 1.0 / 500    #微小時間
const dx = 1.0 / x_max  #微小長さ
const c = 0.5   #伝播速度
const r = c^2 * dt^2 / (2 * dx^2)

function main()
    global t_max
    global x_max
    global dt
    global dx
    global r

    u_init = zeros(x_max + 1)   #初期条件
    u_curr = zeros(x_max + 1)   #U(j)
    u_prev = zeros(x_max + 1)   #U(j-1)
    u_hist = zeros(t_max + 1, x_max + 1)    #記録用

    #Mを定義
    M = zeros(x_max + 1, x_max + 1)
    for x0 = 1:x_max+1
        for x1 = 1:x_max+1
            M[x0, x1] = iseven(x0 + x1) ? 2 * r + 1 : -r
        end
    end
    M = SymTridiagonal(M)

    #初期変位(孤立波)
    w_len = 10#孤立波の幅
    w_pos = 20#孤立波の位置
    for x = 1:w_len+1
        u_init[x+w_pos] = 1 - cos(2 * pi * (x - 1) / w_len)
    end

    for x = 1:x_max+1
        u_curr[x] = u_init[x]
        u_prev[x] = u_init[x]
    end
    @time for t = 1:t_max+1
        u_next = M \ (2 * u_curr - M * u_prev)  #次ステップを計算

        #変数を更新
        for x = 1:x_max+1
            u_hist[t, x] = u_curr[x]
            u_prev[x] = u_curr[x]
            u_curr[x] = u_next[x]
        end
    end

    heatmap(u_hist)
    savefig("waveeq_implicit.png")
    println("end")
end
main()

調べていたら「分数階波動方程式」というのもあるらしく、その計算法もあった。しかし読み解く気力がない。 https://www.researchgate.net/publication/259623599_Numerical_Solution_of_Fractional_Wave_Equation_using_Crank-Nicholson_Method

Karpus-Strong AlgorithmをJuliaで実装

Karpus-Strongアルゴリズム(K-Sアルゴリズム)を↓ページのブロック線図に基づきJuliaで実装した. https://ccrma.stanford.edu/~jos/pasp/Karplus_Strong_Algorithm.html

実行すると生成波形のグラフが表示され音が再生される.

↓生成波形例

f:id:kazmus:20210801133611p:plain
K-Sアルゴリズムによる生成波形の例

以下ソース.

using Random
using Plots
using WAV
using FFTW

mutable struct DelayLine
    _len::Int
    mem::Vector{Float64}
    ptr::Int
    dl_io_do::Function

    DelayLine(_len) = new(_len, zeros(_len), 1)
end

function dl_io_do(dl, mem_in)
    ptr_next = dl.ptr % dl._len + 1
    mem_out = dl.mem[ptr_next]
    dl.mem[dl.ptr] = mem_in
    dl.ptr = ptr_next
    return mem_out
end

function main()
    fs = 8000
    t_snd = 1
    n_delay = 60
    delay1 = DelayLine(n_delay)
    z1 = DelayLine(1)

    rng = MersenneTwister()
    wnoise = randn(rng, Float64, delay1._len)
    delay1.mem = wnoise ./ maximum(abs, wnoise)

    # display(plot(delay1.mem, show = true))

    output = 0
    n_calc = fs * t_snd
    outputs = zeros(n_calc)
    for _i = 1:n_calc
        d1o = dl_io_do(delay1, output)
        z1o = dl_io_do(z1, d1o)
        output = d1o / 2 + z1o / 2
        outputs[_i] = output
    end

    wavplay(outputs, fs)

    display(plot(outputs))
    savefig("KS_wave.png")
end

main()

上記コードの修正

main関数内のコードを以前はトップレベルに書いていた("function main()"とそれに対応する"end"と"main()"の行が無い状態)が、

それを実行すると↓のようなエラーが出たから、関数化してみたら動いた。

この修正方法があっているかわからないが動いたからヨシ

LoadError: UndefVarError: output not defined
in expression starting at D:\Users\USERNAME\Documents\Julia\KS.jl:39
top-level scope at KS.jl:40
eval at boot.jl:360 [inlined]
include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String) at loading.jl:1116

IDEAtomでjunoを導入していたが、この問題に関係あるかは不明。

Python+PyQt5で4次元図形の描画

HSP(Hot Soup Processor)で4次元図形を回転させて描画 - 雑感等

これの二番煎じだけどクラスの設計はMVCを意識したつもり.

4次元図形の回転は上の記事と同様

MVCについて参考になったサイト

動作画面

一番上の行にある3つのスライダーは,4D->3Dに射影するカメラの角度.

その下の行にある2つのスライダーは,3D->2Dに射影するカメラの角度.

その下には描画された図形.この描画範囲内でホイールをコロコロすると,図形が拡大縮小する.

f:id:kazmus:20210526102212p:plain

pythonソース

from PyQt5.QtWidgets import QWidget, QApplication, QSlider, QGraphicsView, QGraphicsScene, QVBoxLayout, QHBoxLayout, \
    QGraphicsItem
from PyQt5.QtGui import QPen
from PyQt5.QtCore import Qt, QRectF
import sys
import numpy as np


class Model(object):
    def __init__(self):
        self.lines = None
        self._init_lines()
        pass

    def _init_lines(self):
        """
        線(点座標の対)を初期化
        点座標は4次元
        :return:
        """
        self.angle4d = [0, 0, 0]
        self.angle3d = [0, 0]
        p_o = np.array([0, 0, 0, 0])
        p_x = np.array([1, 0, 0, 0])
        p_y = np.array([0, 1, 0, 0])
        p_z = np.array([0, 0, 1, 0])
        p_w = np.array([0, 0, 0, 1])
        p_any = np.array([1, 1, 1, 1])
        p_any = p_any / np.linalg.norm(p_any)
        self.lines = np.stack([
            np.vstack([p_o, p_x]),
            np.vstack([p_o, p_y]),
            np.vstack([p_o, p_z]),
            np.vstack([p_o, p_w]),
            np.vstack([p_o, p_any])
        ])
        pass

    def set_t(self, val):
        self.angle4d[0] = val
        pass

    def set_u(self, val):
        self.angle4d[1] = val
        pass

    def set_v(self, val):
        self.angle4d[2] = val
        pass

    def set_t3(self, val):
        self.angle3d[0] = val
        pass

    def set_u3(self, val):
        self.angle3d[1] = val
        pass

    def _calc_projection4d(self, point_4, angle4d_cam, angle3d_cam):
        """
        投影を計算
        :param point_4:
        :param angle:
        :return:
        """
        t, u, v = angle4d_cam
        t3, u3 = angle3d_cam

        st = np.sin(t)
        ct = np.cos(t)
        su = np.sin(u)
        cu = np.cos(u)
        sv = np.sin(v)
        cv = np.cos(v)

        st3 = np.sin(t3)
        ct3 = np.cos(t3)
        su3 = np.sin(u3)
        cu3 = np.cos(u3)

        proj_4to3 = np.array([
            [-st, ct, 0, 0],
            [-ct * su, -st * su, cu, 0],
            [-ct * cu * sv, -st * cu * sv, -su * sv, cv]
        ])
        proj_3to2 = np.array([
            [-st3, ct3, 0],
            [-ct3 * su3, -st3 * su3, cu3],
        ])
        return proj_3to2 @ proj_4to3 @ point_4
        pass

    def projection(self):
        """
        投影された座標データを返す
        :return:
        """
        angle4d = np.array(self.angle4d)
        angle3d = np.array(self.angle3d)
        # ↓でやりたいことは,行列の第2方向(第0,第1に次ぐ第2方向)の各要素(4要素のベクトル)に対して,_calc_projection4dを適用したい
        return np.apply_along_axis(lambda x: self._calc_projection4d(x, angle4d, angle3d), 2, self.lines)
        pass

    pass


class Plotter(QGraphicsItem):
    def __init__(self, width, height, model):
        super(Plotter, self).__init__()
        self.width = width
        self.height = height
        self.model = model
        self.zoom = 200

        pass

    def plottable_centre(self):
        return np.array([self.width / 2, self.height / 2, self.width / 2, self.height / 2])
        pass

    def do_plot(self):
        lines = self.model.projection()
        centre = self.plottable_centre()
        self.plottable_lines = list(map(lambda x: list(map(int, x.flatten() * self.zoom + centre)), lines))
        self.update()
        pass

    def paint(self, painter, option, widget):
        colours = [Qt.red, Qt.green, Qt.blue, Qt.black, Qt.darkCyan, Qt.darkYellow, Qt.darkYellow]
        line_colour = zip(self.plottable_lines, colours)
        for line, colour in line_colour:
            pen = QPen()
            pen.setWidth(3)
            pen.setBrush(colour)
            painter.setPen(pen)
            painter.drawLine(*line)
            pass
        pass

    def boundingRect(self):
        return QRectF(0, 0, self.width, self.height)
        pass

    def wheelEvent(self, event):
        d = event.delta()
        self.zoom += d / 10
        self.do_plot()
        pass

    pass


class View(QWidget):
    def __init__(self, model):
        super(View, self).__init__()
        self.model = model
        pass

    def register(self, controller):
        self.controller = controller
        self.init_ui()
        pass

    def init_ui(self):
        """
        ウィジェットの配置
        :return:
        """
        self.setWindowTitle("4D rotation")

        self.scene_size = 500
        self.plotter = Plotter(self.scene_size, self.scene_size, self.model)

        controllers_4d_layout = QHBoxLayout()

        controllers_4d = self.controller.controllers_4d
        for ix in range(len(controllers_4d)):
            sld = QSlider(Qt.Horizontal, self)
            sld.setValue(33)
            controllers_4d_layout.addWidget(sld)
            sld.valueChanged.connect(controllers_4d[ix])
            controllers_4d[ix](sld.value())
            pass

        controllers_3d_layout = QHBoxLayout()
        controllers_3d = self.controller.controllers_3d
        for ix in range(len(controllers_3d)):
            sld = QSlider(Qt.Horizontal, self)
            sld.setValue(33)
            controllers_3d_layout.addWidget(sld)
            sld.valueChanged.connect(controllers_3d[ix])
            controllers_3d[ix](sld.value())
            pass

        self.graphics_view = QGraphicsView()
        scene = QGraphicsScene(self.graphics_view)
        scene.setSceneRect(0, 0, self.scene_size, self.scene_size)
        self.graphics_view.setScene(scene)
        scene.addItem(self.plotter)

        # self.setLayout(grid)
        main_layout = QVBoxLayout()
        main_layout.setAlignment(Qt.AlignTop)
        main_layout.addLayout(controllers_4d_layout)
        main_layout.addLayout(controllers_3d_layout)
        main_layout.addWidget(self.graphics_view)
        self.setLayout(main_layout)

        pass

    pass


class Controller(object):
    def __init__(self, view, model):
        self.view = view
        self.model = model

        self.controllers_4d = [self.chg_sld_t, self.chg_sld_u, self.chg_sld_v]
        self.controllers_3d = [self.chg_sld_t3, self.chg_sld_u3]
        self.view.register(self)
        self.do_plot()
        pass

    def do_plot(self):
        self.view.plotter.do_plot()
        pass

    def chg_sld_t(self, value):
        self.model.set_t(value * 2 * np.pi / 100)
        self.do_plot()
        pass

    def chg_sld_u(self, value):
        self.model.set_u(value * 2 * np.pi / 100)
        self.do_plot()
        pass

    def chg_sld_v(self, value):
        self.model.set_v(value * 2 * np.pi / 100)
        self.do_plot()
        pass

    def chg_sld_t3(self, value):
        self.model.set_t3(value * 2 * np.pi / 100)
        self.do_plot()
        pass

    def chg_sld_u3(self, value):
        self.model.set_u3(value * 2 * np.pi / 100)
        self.do_plot()
        pass

    pass


def main():
    app = QApplication(sys.argv)
    model = Model()
    view = View(model)
    controller = Controller(view, model)

    view.show()
    sys.exit(app.exec_())
    pass


if __name__ == '__main__':
    main()
    pass

PC変えてからHSPインタプリタ入れてなくて,入れるのめんどくさかった.

pythonでのgui開発とオブジェクト指向っぽい書き方の勉強がてらpythonに移植?した.

PlotterクラスとViewクラスに分けて実装した(参考にしたサイト・ソース↓がQGraphicsItemのクラスとViewクラスに分けてた)がMVC的には良かったのかわからない.

pythonでmidiメッセージ送信

pythonのライブラリmidoを使用してリアルタイムにmidiメッセージを送信する.

以下のプログラムでは,ランダムな音高・音価・サステインペダルの踏み度合いを送信する.

pianoteqのstandaloneとloopMIDIも合わせて使っている.

import mido
import time
import numpy as np

ports = mido.get_output_names()
midi_port_name = ports[1]
time.sleep(1)
with mido.open_output(midi_port_name) as outport:
    outport.reset()
    time.sleep(1)


    def sustain_off():
        outport.send(mido.Message("control_change", channel=0, control=64, value=0, time=1))
        pass


    def sustain_on():
        outport.send(mido.Message("control_change", channel=0, control=64, value=127, time=1))
        pass


    def sustain(value):
        outport.send(mido.Message("control_change", channel=0, control=64, value=value, time=1))
        pass


    sustain_on()
    for n in range(1000):
        print(n)
        note = np.random.randint(21, 108 + 1)
        # vel_a = 2  # 夜用
        # vel_b = 12  # 夜用
        vel_a = 1.5
        vel_b = 1.5
        velo_on = round(np.random.beta(vel_a, vel_b) * 127)  # 朝用
        velo_off = round(np.random.beta(vel_a, vel_b) * 127)  # 朝用
        # dura = np.random.gamma(16 / 1.5, 0.2 * 1.5)  # 夜用
        dura = np.random.gamma(1, 0.2)  # 朝用
        sust = round(np.random.beta(0.8, 0.2) * 127)

        sustain(sust)
        time.sleep(0.01)
        outport.send(mido.Message("note_on", channel=0, note=note, velocity=velo_on, time=1))
        time.sleep(dura)
        outport.send(mido.Message("note_off", channel=0, note=note, velocity=velo_off, time=1))
        pass
    time.sleep(3)
    sustain_off()

    time.sleep(1)
    outport.reset()
    pass