扇型に内接する楕円を求めた
――かったが、実質的には二等辺三角形に内接する楕円の式が求まった。
扇型は、半径を1、中心角をとする。
楕円は、一方(扇型の径方向)の半径を、もう一方の半径をとする。
各図形は図の座標に配置し、とを既知の値とすると、は次式で与えられる。
このを使い、上図の楕円の式は次式で与えられる。
このページで式を動かせる:二等辺三角形に内接する楕円 – GeoGebra
扇型に内接する楕円を求めた
――かったが、実質的には二等辺三角形に内接する楕円の式が求まった。
扇型は、半径を1、中心角をとする。
楕円は、一方(扇型の径方向)の半径を、もう一方の半径をとする。
各図形は図の座標に配置し、とを既知の値とすると、は次式で与えられる。
このを使い、上図の楕円の式は次式で与えられる。
このページで式を動かせる:二等辺三角形に内接する楕円 – GeoGebra
[^\(\)]+(?=\([^\(\)]+\))|[^\(\)]+$|^[^\(\)]+
この正規表現は、orを意味するパイプ記号"|"と、3つの部分(下記ABC)で構成されている。
A | B | C
↓
[^\(\)]+(?=\([^\(\)]+\)) | [^\(\)]+$ | ^[^\(\)]+
つまり、「AまたはBまたはC」のパターンを表す。
これは、さらに3つの要素に分けられる。
D E F
↓
[^\(\)] + $
最初の角括弧[]は、「指定した文字を含まない」という意味の表現:[^~~~~]だ。
"[^"と"]"に挟まれた各文字以外の文字にマッチする。ここでいう各文字は \( と \)だ。エスケープされているが要は丸括弧"("または")"を意味する。
つまり、[^\(\)]は「丸括弧でない何らかの1文字」を意味する。
次、+は「前の文字が1文字以上連続する」を意味する。
次、$は、「行末」を意味する。
ここまで合わせて考えると[^\(\)]+$は「丸括弧以外の文字が1文字以上連続し、その直後に行末が来る」パターンを表す。
簡単に言うと、「丸括弧から行末までの文字 ※丸括弧は含まず」にマッチする。
例えば入力が「abc0123(hijk」なら、「hijk」にマッチする。
最初の^は、「行頭」を表す。※"[^"の"^"とは意味が異なる。
次は上で見た、[^\(\)]+と全く同じだ。
つまり、^[^\(\)]+は「行頭の直後に、丸括弧以外の文字が1文字以上連続する」パターンを表す。
簡単に言うと、「行頭から丸括弧までの文字 ※丸括弧は含まず」にマッチする。
例えば入力が「abc0123(hijk」なら、「abc0123」にマッチする。
まず、大枠は
G (?= H )
↓
[^\(\)]+ (?= \([^\(\)]+\) )
G (?= H ) のパターンは、肯定先読みと呼ばれる。
意味は、「G。ただし、Hが直後にあるもの」のパターンを表す。
「マッチしてほしいのはGだけど、Gならなんでもいいわけじゃなくて、Hの直前にあるGだけにマッチしてほしい」という事。
さて、GとHを見ていく。
Gは上で見た[^\(\)]+と全く同じだ。 つまり、「丸括弧以外の文字が1文字以上連続する」パターンを表す。
次、H: \([^\(\)]+\)は、さらに3つの要素に分解できて、 \( [^\(\)]+ \)となる。
\(は開き括弧、 [^\(\)]+は「丸括弧以外の文字が1文字以上連続する」パターン、 \)は閉じ括弧。
Hはつまり、「丸括弧で囲まれた文字列(ただし囲んでいる丸括弧も含む)」パターンを表す。
つまり[^\(\)]+ (?= \([^\(\)]+\) )は、 「丸括弧で囲まれた文字列の直前にある、丸括弧以外の文字列」のパターンを表す。
例えば入力が「)abcd(efgh)0123(4567)ijkl」なら、「abcd」と「0123」にマッチする。
いろいろ調整してみたら自励振動させられてないっぽい。管体部分でループ(ハウリングみたいな状態)していただけだった。
物理モデリング音源の特許(だいたい権利失効してる)の実施例にある、励振部をpythonで実装した。
一応音は出るが、パラメータの調整が甘くか弱い音しか鳴らない。
とりあえず音が出る状態でメモしておく。
■pythonソース
import pyaudio import numpy as np from collections import deque import matplotlib.pyplot as plt import scipy.signal as signal fs_sound = 44100 # [1/s] sampling freq. fs_internal = 44100 c = 340 # [m/s] speed of sound class DelayLine: def __init__(self, n_delay): n_delay = np.round(n_delay).astype(np.int64) self._dl = deque([0] * n_delay) pass def update(self, x_in): self._dl.append(x_in) x_out = self._dl.popleft() return x_out pass class Filter(DelayLine): def __init__(self, coef_fir): n_taps = coef_fir.size super().__init__(n_taps) self.coef_fir = coef_fir[::-1] pass def update(self, x_in): self._dl.append(x_in) self._dl.popleft() x_out = self._dl @ self.coef_fir return x_out pass def sound_play(samples): p = pyaudio.PyAudio() stream = p.open(format=pyaudio.paFloat32, channels=1, rate=fs_sound, frames_per_buffer=1024, output=True) t_end = 1 n_samples = t_end * fs_sound n_t = np.arange(0, fs_sound) f = 440 if np.max(np.abs(samples)) > 1: samples = samples / np.max(np.abs(samples)) * 0.9 pass stream.write(samples.astype(np.float32).tobytes()) stream.close() pass def filter_analy(a, fs): freq, res = signal.freqz(a) tate = 3 yoko = 1 plt.subplot(tate, yoko, 1) plt.plot(0.5 * fs * freq / np.pi, np.abs(res)) plt.subplot(tate, yoko, 2) plt.plot(0.5 * fs * freq / np.pi, np.angle(res)) plt.subplot(tate, yoko, 3) plt.plot(a[::-1], ".-") plt.show() pass def main(): # FILTER SETTING ########################################### fs = fs_internal fc = 2000 freqs = np.array([0, fc, fc + 1000, fs / 2]) a = signal.firwin2(numtaps=500, freq=freqs / (fs / 2), gain=[1, 1, 0, 0], nfreqs=fs) a = signal.minimum_phase(a, "homomorphic") fl1 = Filter(a) fl_press2reed = Filter(a) fc = 30 freqs = np.array([0, fc, fc + 50, fs / 2]) a = signal.firwin2(numtaps=2001, freq=freqs / (fs / 2), gain=[0, 0, 1, 1], nfreqs=fs) a = signal.minimum_phase(a, "homomorphic") fl_outhpf = Filter(a) ############################################ dl_len = 0.8 # [m] delay line length n_dl = fs_internal * dl_len / c print(n_dl) # filter_analy(a, fs) dl1 = DelayLine(n_dl) dl2 = DelayLine(n_dl) dl_s = DelayLine(1) t_simulation = 2 n_sim = fs_internal * t_simulation xs_input = np.zeros(n_sim) in_max = 3 in_min = 0.15 n_on1 = 10000 n_on2 = 20000 n_off1 = 30000 n_off2 = 80000 # xs_input[0] = 1 for n in range(n_sim): if n_on1 <= n and n < n_on2: x = in_min + (in_max - in_min) / (n_on2 - n_on1) * (n - n_on1) pass elif n_on2 <= n and n < n_off1: x = in_max pass elif n_off1 <= n and n < n_off2: x = in_min + (in_max - in_min) * (1 - 1 / (n_off2 - n_off1) * (n - n_off1)) pass else: x = in_min pass xs_input[n] = x pass xs_output = np.zeros(n_sim) xs_reed = np.zeros(n_sim) xs_volvel = np.zeros(n_sim) k_refl_reed = 0.1 # INIT ############################ do1 = 0.0 do2 = 0.0 do_s = 0.0 def press2reed(x): out = x if x < -1: out = -1 pass elif x > 1: out = 1 pass else: out = x pass # return 1 / (1 + np.exp(-x)) return out pass def press2vel(x): # return np.abs(x) return x pass def volvel2press(x): return x pass for n in range(n_sim): x_n = xs_input[n] s_press = do_s - x_n s_press_reed = s_press - 0.8 s_reed_area = press2reed(fl_press2reed.update(s_press_reed)) # s_reed_area = press2reed(s_press_reed) xs_reed[n] = s_reed_area s_velocity = press2vel(s_press) s_volvel = s_reed_area * s_velocity s_reedout = volvel2press(s_volvel) xs_volvel[n] = s_volvel s_forward_d = (1 + k_refl_reed) * s_reedout + (-k_refl_reed) * do2 do1 = dl1.update(s_forward_d) x_n = fl_outhpf.update(do1) s_backward_1 = fl1.update(x_n) do2 = dl2.update(s_backward_1) s_backward_2 = (1 - k_refl_reed) * do2 + k_refl_reed * s_reedout do_s = dl_s.update(s_backward_2) xs_output[n] = x_n pass plt.plot(xs_output, label="out") plt.plot(xs_input, "--", label="in") plt.plot(xs_reed, label="reed") plt.plot(xs_volvel, label="volvel") plt.legend() plt.show() sound_play(xs_output) pass if __name__ == '__main__': main()
出典
#
配列の要素数はコンパイル時に決定されているため、入力に応じて長さを伸縮できない。
下記はコンパイルエラーになる。(要素数arr_lenが変数のためコンパイルエラー)
let arr_len: usize = 8; let arr = [0 as u8; arr_len];
一方、Vecだとコンパイル時に要素数が決定されていなくても良い。
下記のようにvec!マクロで宣言すると楽。(要素数はvec_len)
let vec_array = vec![0 as u8; vec_len];
#
下記rustソースでは、
標準入力に応じてvec_arrayの要素数が変化するのが確認でき、 実行時に(入力に応じて)要素数を決定できる。
use std::io; fn main() { //文字列格納先 let mut input_str = String::new(); //標準入力から文字列を読む println!("array len:"); io::stdin() .read_line(&mut input_str) .expect("failed to read line"); //usize型に変換 let vec_len: usize = input_str.trim().parse().expect("failed to parse"); //u8型を要素に持つVec配列を生成 let vec_array = vec![0 as u8; vec_len]; //生成した配列長と、配列を表示 println!("array len={}\tarray={:?}", vec_array.len(), vec_array); }
実行例(10を入力)
array len: 10 array len=10 array=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0] Process finished with exit code 0
C言語で動的配列と調べるとmallocがすぐに出るが、 rustだと検索してもvec!が直球で出てこなかった気がした。(ドキュメントをちゃんと読めば書いてあると思うが)
↓このあたりのソースを改変して、圧縮後のzipファイルの保存先を指定できるようにした。
batプログラムソースは以下。
【保存先フォルダ(フルパス)】のところに、保存先フォルダのパスを書く。
ZIP_PATHには7zipの実行ファイルのパスを指定する。
@echo off set ZIP_PATH="C:\Program Files\7-Zip\7z.exe" for %%f in (%*) do ( call :filename %%f ) exit /b :filename %ZIP_PATH% a -tzip "【保存先フォルダ(フルパス)】\%~n1.zip" "%~1"
%%fから、圧縮するフォルダ名だけ切り出す方法がわからなかったので、
サブルーチン:filenameを呼び出して、%~n1でフォルダ名だけ切り出すようにした。
バグを直した&メッシュ数に応じて出力画像の比率を変えるようにした。
■出力画像
■Rustソース
use std::time::{Duration, Instant}; use plotters::prelude::*; const PLT_RATIO: usize = 20; const T_MAX: usize = 1000; const X_MAX: usize = 10; const Y_MAX: usize = 31; const M_ELEM: usize = 3 * 4 + 4 * 2 * (X_MAX - 2) + 4 * 2 * (Y_MAX - 2) + 5 * (X_MAX - 2) * (Y_MAX - 2); const T_WINDOW: usize = 4; struct SparseMatrix { col: [usize; M_ELEM], row: [usize; M_ELEM], val: [f64; M_ELEM], } fn sub_main() { println!("M_ELEM = {:?}", M_ELEM); // plot //////////////////////////////////////////////////////////////////////////////////////// let root = BitMapBackend::gif("testrustimage.gif", ((X_MAX * PLT_RATIO) as u32, (Y_MAX * PLT_RATIO) as u32), 50).unwrap() .into_drawing_area(); let cells = root.split_evenly((Y_MAX, X_MAX)); let mut cell_nx; let mut cell_ny; let mut val; let mut col; // calc //////////////////////////////////////////////////////////////////////////////////////// const DT: f64 = 1.0 / 44100.0; const DX: f64 = 0.05; const C: f64 = 347.0; const COURANT: f64 = C * DT / DX; println!("DT:{},DX:{},C:{},COURANT:{}", DT, DX, C, COURANT); //u[time][y][x] let mut u: [[[f64; X_MAX]; Y_MAX]; T_WINDOW] = [[[0.0f64; X_MAX]; Y_MAX]; T_WINDOW]; let mut M: [[f64; X_MAX * Y_MAX]; X_MAX * Y_MAX] = [[0.0f64; X_MAX * Y_MAX]; X_MAX * Y_MAX]; let mut x: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX]; let mut b: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX]; let r = 0.5; let nu2 = COURANT * COURANT; //define M { let diag = 4.0 * r * nu2 + 2.0; let oth_elem = -r * nu2; let mut curr_ind = 0; for ny in 0..Y_MAX { for nx in 0..X_MAX { curr_ind = ny * X_MAX + nx; M[curr_ind][curr_ind] = diag; if nx > 0 { M[curr_ind][curr_ind - 1] = oth_elem; } if nx < X_MAX - 1 { M[curr_ind][curr_ind + 1] = oth_elem; } if ny > 0 { M[curr_ind][curr_ind - X_MAX] = oth_elem; } if ny < Y_MAX - 1 { M[curr_ind][curr_ind + X_MAX] = oth_elem; } } } } let M = M; println!("u:{}", u[0][0][X_MAX - 1]); let mut u_max = 2.0; let mut u_max_t: usize = 0; let mut u_max_x: usize = 0; let mut u_max_y: usize = 0; let mut u_min = -1.0; let mut u_min_t: usize = 0; let mut u_min_x: usize = 0; let mut u_min_y: usize = 0; for nx in 1..X_MAX - 1 { u[1][2][nx] = 1.0; } let mut t_curr; let mut t_prev; let mut t_next; let loop_start = Instant::now(); //main loop //gauss seidel: solve Mx=b { let omega = 1.0; let mut x_temp = 0.0; let mut m_temp = 0.0; let mut gs_diff = 0.0; let mut curr_ind = 0; for nt in 1..T_MAX - 1 { t_curr = nt % T_WINDOW; t_prev = (nt - 1) % T_WINDOW; t_next = (nt + 1) % T_WINDOW; //define b for ny in 1..=Y_MAX - 2 { for nx in 1..=X_MAX - 2 { curr_ind = ny * X_MAX + nx; b[curr_ind] = r * nu2 * (u[t_prev][ny][nx - 1] + u[t_prev][ny - 1][nx] + u[t_prev][ny][nx + 1] + u[t_prev][ny + 1][nx]) + 2.0 * (1.0 - r) * nu2 * (u[t_curr][ny][nx - 1] + u[t_curr][ny - 1][nx] + u[t_curr][ny][nx + 1] + u[t_curr][ny + 1][nx]) - (4.0 * r * nu2 + 2.0) * u[t_prev][ny][nx] + (8.0 * (r - 1.0) * nu2 + 4.0) * u[t_curr][ny][nx]; } } //gauss seidel for lp in 0_usize..200_usize { gs_diff = 0.0; for ny in 0..X_MAX * Y_MAX { for nx in 0..ny { m_temp -= M[ny][nx] * x[nx]; } for nx in (ny + 1)..X_MAX * Y_MAX { m_temp -= M[ny][nx] * x[nx]; } x_temp = (b[ny] + m_temp) / M[ny][ny]; gs_diff += (x_temp - x[ny]).abs(); // println!("{:?}:{:?}:{:?}", nt, lp, x_temp - x[ny]); x[ny] = omega * x_temp + (1.0 - omega) * x[ny]; m_temp = 0.0; } if gs_diff < f64::MIN_POSITIVE { println!("break {}:{}", nt, lp); break; } } println!("{:?}:{:?}", nt, gs_diff); //x->u for ny in 0..Y_MAX { for nx in 0..X_MAX { curr_ind = ny * X_MAX + nx; u[t_next][ny][nx] = x[curr_ind]; } } // plot //////////////////////////////////////////////////////////////////////////// for (area, num) in (&cells).into_iter().zip(0..) { cell_nx = num % X_MAX; cell_ny = num / X_MAX; val = u[t_curr][cell_ny][cell_nx]; col = (val - u_min) / (u_max - u_min); area.fill(&RGBColor((col * 255.0) as u8, (col * 255.0) as u8, (col * 255.0) as u8)).unwrap(); } root.present().expect("unable to write"); } } let duration = loop_start.elapsed(); println!("elapsed:{:?}", duration); println!("max u[{}][{}][{}]={}", u_max_t, u_max_y, u_max_x, u_max); println!("min u[{}][{}][{}]={}", u_min_t, u_min_y, u_min_x, u_min); } fn main() { { const STACK_SIZE: usize = 1024 * 1024 * 1024 * 1024; std::thread::Builder::new() .stack_size(STACK_SIZE) .spawn(sub_main) .unwrap() .join() .unwrap(); } }
※バグあり メッシュ数を縦≠横にするとおかしくなる
とりあえずソースと結果を貼っておく。
陰解法と書いたが、実際は陰解法と陽解法を平均した差分方程式を使っている。(クランクニコルソン法っぽくした)
↓の記事で導出した式と、陽解法を組み合わせている。
波動方程式の有限差分法を半自動で導出(sympy) - 雑感等
ソースの条件だとクーラン数=0.1573<1で陽解法でも解ける。
const DT: f64 = 1.0 / 44100.0;
const DX: f64 = 0.05;
このあたりを変えるとクーラン数が変わる。
クーラン数を1以上にしても一応解けたが、上げるとガウス・ザイデル法の反復回数が増える(収束が遅くなる)。
下のソースだと反復回数に上限を設けているから、収束し無くても打ち切る場合がある。
ある程度まで反復回数が増えて、一定を超えると、誤差が増えていくはず。
■出力結果
■Rustプログラムのソース
use std::time::{Duration, Instant}; use plotters::prelude::*; const T_MAX: usize = 100; const X_MAX: usize = 50; const Y_MAX: usize = 50; const M_ELEM: usize = 3 * 4 + 4 * 2 * (X_MAX - 2) + 4 * 2 * (Y_MAX - 2) + 5 * (X_MAX - 2) * (Y_MAX - 2); const T_WINDOW: usize = 4; struct SparseMatrix { col: [usize; M_ELEM], row: [usize; M_ELEM], val: [f64; M_ELEM], } fn sub_main() { println!("M_ELEM = {:?}", M_ELEM); // plot //////////////////////////////////////////////////////////////////////////////////////// let root = BitMapBackend::gif("testrustimage.gif", (200, 200), 50).unwrap() .into_drawing_area(); let cells = root.split_evenly((X_MAX, Y_MAX)); let mut cell_nx; let mut cell_ny; let mut val; let mut col; // calc //////////////////////////////////////////////////////////////////////////////////////// const DT: f64 = 1.0 / 44100.0; const DX: f64 = 0.05; const C: f64 = 347.0; const COURANT: f64 = C * DT / DX; println!("DT:{},DX:{},C:{},COURANT:{}", DT, DX, C, COURANT); //u[time][y][x] let mut u: [[[f64; X_MAX]; Y_MAX]; T_WINDOW] = [[[0.0f64; X_MAX]; Y_MAX]; T_WINDOW]; let mut M: [[f64; X_MAX * Y_MAX]; X_MAX * Y_MAX] = [[0.0f64; X_MAX * Y_MAX]; X_MAX * Y_MAX]; let mut x: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX]; let mut b: [f64; X_MAX * Y_MAX] = [0.0f64; X_MAX * Y_MAX]; let r = 0.5; let nu2 = COURANT * COURANT; //define M { let diag = 4.0 * r * nu2 + 2.0; let oth_elem = -r * nu2; let mut curr_ind = 0; for ny in 0..Y_MAX { for nx in 0..X_MAX { curr_ind = ny * X_MAX + nx; M[curr_ind][curr_ind] = diag; if nx > 0 { M[curr_ind][curr_ind - 1] = oth_elem; } if nx < X_MAX - 1 { M[curr_ind][curr_ind + 1] = oth_elem; } if ny > 0 { M[curr_ind][curr_ind - X_MAX] = oth_elem; } if ny < Y_MAX - 1 { M[curr_ind][curr_ind + X_MAX] = oth_elem; } } } } let M = M; println!("u:{}", u[0][0][X_MAX - 1]); let mut u_max = 2.0; let mut u_max_t: usize = 0; let mut u_max_x: usize = 0; let mut u_max_y: usize = 0; let mut u_min = -1.0; let mut u_min_t: usize = 0; let mut u_min_x: usize = 0; let mut u_min_y: usize = 0; u[1][10][20] = 1.0; let mut t_curr; let mut t_prev; let mut t_next; let loop_start = Instant::now(); //main loop //gauss seidel: solve Mx=b { let omega = 1.0; let mut x_temp = 0.0; let mut m_temp = 0.0; let mut gs_diff = 0.0; let mut curr_ind = 0; for nt in 1..T_MAX - 1 { t_curr = nt % T_WINDOW; t_prev = (nt - 1) % T_WINDOW; t_next = (nt + 1) % T_WINDOW; //define b for ny in 1..=Y_MAX - 2 { for nx in 1..=X_MAX - 2 { curr_ind = ny * X_MAX + nx; b[curr_ind] = r * nu2 * (u[t_prev][ny][nx - 1] + u[t_prev][ny - 1][nx] + u[t_prev][ny][nx + 1] + u[t_prev][ny + 1][nx]) + 2.0 * (1.0 - r) * nu2 * (u[t_curr][ny][nx - 1] + u[t_curr][ny - 1][nx] + u[t_curr][ny][nx + 1] + u[t_curr][ny + 1][nx]) - (4.0 * r * nu2 + 2.0) * u[t_prev][ny][nx] + (8.0 * (r - 1.0) * nu2 + 4.0) * u[t_curr][ny][nx]; } } //gauss seidel for lp in 0_usize..200_usize { gs_diff = 0.0; for ny in 0..X_MAX * Y_MAX { for nx in 0..ny { m_temp -= M[ny][nx] * x[nx]; } for nx in (ny + 1)..X_MAX * Y_MAX { m_temp -= M[ny][nx] * x[nx]; } x_temp = (b[ny] + m_temp) / M[ny][ny]; gs_diff += (x_temp - x[ny]).abs(); // println!("{:?}:{:?}:{:?}", nt, lp, x_temp - x[ny]); x[ny] = omega * x_temp + (1.0 - omega) * x[ny]; m_temp = 0.0; } if gs_diff < f64::MIN_POSITIVE { println!("break {}:{}", nt, lp); break; } } println!("{:?}:{:?}", nt, gs_diff); //x->u for ny in 0..Y_MAX { for nx in 0..X_MAX { curr_ind = ny * X_MAX + nx; u[t_next][ny][nx] = x[curr_ind]; } } // plot //////////////////////////////////////////////////////////////////////////// for (area, num) in (&cells).into_iter().zip(0..) { cell_nx = num % X_MAX; cell_ny = num / X_MAX; val = u[t_curr][cell_nx][cell_ny]; col = (val - u_min) / (u_max - u_min); area.fill(&RGBColor((col * 255.0) as u8, (col * 255.0) as u8, (col * 255.0) as u8)).unwrap(); } root.present().expect("unable to write"); } } let duration = loop_start.elapsed(); println!("elapsed:{:?}", duration); println!("max u[{}][{}][{}]={}", u_max_t, u_max_y, u_max_x, u_max); println!("min u[{}][{}][{}]={}", u_min_t, u_min_y, u_min_x, u_min); } fn main() { { const STACK_SIZE: usize = 1024 * 1024 * 1024 * 1024; std::thread::Builder::new() .stack_size(STACK_SIZE) .spawn(sub_main) .unwrap() .join() .unwrap(); } }