Skip to content

phasetr/mathcodes

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

mathcodes

I collected my scientific codes, especially math and physics.

雑多な覚え書き

コードへのコメント

Julia > samples > movie_0001.jl など, いくつかのファイルはこの readme に YouTube の対応動画へのリンクを貼っている. これは途中からはじめたのだが, わざわざ分離して readme に書くのも面倒になったので, ファイルに直接 URL を書くのを許すことにした. この目的なら gif にして Jupyter notebook を使うのも一手と思うが, Jupyter notebook は json で容量があり, 差分チェックもはてしなく鬱陶しい点で, 2020-04 時点では回避している.

どんなときに何をどうするといいかはいろいろ検討している.

python

Jupyter がすごいのでそれを使ったコードも書いてみたい. ここに行けばネット越しに Jupyter が使えるので, インストールしたりコマンド叩いたりがつらい人はとりあえずこちらで実験してみよう.

以下, インストールについての注意を書いておく. 2016-08 時点ではQiita のこの記事が参考になる.

詳しくはそちらを見てもらうことにして, ここでは簡単にまとめておく.

readme を書いた時点では次のバージョンで動かしている.

  • Python 3.5.1 :: Anaconda 4.1.0 (x86_64)

ライブラリのバージョンも書いた方がいいのかもしれないが, それはサボる.

Anaconda でのインストール

Windows (未検証)

Mac (検証)

Homebrew で次のコマンド実行でインストールした. brew upgrade でかなり時間がかかった. 注意してほしい.

また次のコマンドリストの # が入っている箇所では, # 以下は入力しないようにすること. コメントのつもりで書いている.

$ brew update && brew upgrade
$ brew install pyenv
$ echo 'export PATH="$HOME/.pyenv/bin:$PATH"' >> ~/.bash_profile
$ echo 'eval "$(pyenv init -)"' >> ~/.bash_profile
$ exec $SHELL -l
$ pyenv install -l | grep anaconda # Check the newest version
$ pyenv install anaconda3-4.1.0    # Input the newest
$ pyenv global anaconda3-4.1.0     # Input the newest
$ python --version

TODO VPython インストール

2016-08 時点で Python3 対応ができていないようだ. 計算物理の文献で VPython を使っているケースを見かけるので, Python3 に対応次第, こちらも追記したい.

Jupyter 起動

シンプルに

本来はコマンドプロンプトやターミナル (いわゆる「黒いやつ」) で起動しないといけない. しかしこのハードルは高いだろうという気もする. そこでダブルクリックすれば自動で起動できるファイルを作った.

  • jupyter/runjupyter.bat
  • jupyter/runjupyter.command

Windows の方は bat ファイルを, Mac の方は command ファイルを使うこと. これを適当なフォルダに置いてダブルクリックで実行すれば, そのフォルダで jupyter が起動する.

詳細

jupyter/fundamental.ipynb にも書いておいたが, そもそも起動できないのに確認しようもない. 簡単に書いておく.

$ cd some_directory # 適当なディレクトリに移動
$ jupyter notebook

これでそのディレクトリ上で jupyter が起動する. ファイルを作ったりすると some_directory にファイルが作られる.

Jupyter Notebook については jupyter/fundamental.ipynb にまとめてあるのでそちら参照. 次のようなことがまとめてある.

  • Markdown 記法が使えること
  • mathjax が使えるので TeX が使えること.
  • 簡単なキーバインドまとめ.
  • Python の基礎の基礎.
  • IPython の基礎の基礎.

簡単なものでもアニメーションやグラフを描けるのが楽しい. jupyter/fundamental.ipynb を軽く眺めたら, 次に jupyter/math_simple_graph.ipynb を見てほしい. あとは適当に書き進めていく.

julia

速くて書くのも楽というので試してみる.

Emacs Julia-Repl

key action
C-c C-a activate if there is a Project.toml in parent directories
C-u C-c C-a activate home project
C-c C-b send whole buffer to REPL (using include)
C-u C-c C-b send whole buffer to REPL (directly)
C-c C-c send region (when applicable) or line to REPL
C-c C-d invoke @doc on symbol
C-c C-e invoke @edit on region (when applicable) or line
C-c C-l list methods of a function
C-c C-m expand macro
C-c C-p change directory to that of the buffer
C-c C-s prompt for buffer name suffix
C-c C-t send whole buffer to REPL (using Revise.includet)
C-c C-v prompt for Julia executable
C-c C-z raise the REPL or create a new one
C-RET send line to REPL (without bracketed paste)

資料メモ

自分用チュートリアル・クックブック用メモ

Plots

メモ

  • Julia では高速に計算したいコードは函数の中に入れること
    • 関数実行時に与えられた引数の型情報を使ってその函数をネイティブコードにコンパイルしてから実行する仕組みがあるため.
    • 最初の実行時にはコンパイルにも時間が取られるので注意すること.
  • Unicode の数学シンボルがタブで打てる
    • Julia の REPL やいくつかの Julia 編集環境で使える: 少なくとも Emacs の julia-mode では使えた
    • 例: δ -> \delta-tab
    • 例: α̂₂ -> \alpha-tab-\hat-tab-_2-tab

samples

movie_0001.jl

plots_tutorial_0002.jl

movie_0003.jl

movie_0006.jl

visualization

VTK

Gnuplot

生成した数値計算結果ファイルの削除

次のコマンドで良さそう.

git clean -dnX ./   # まずこれでチェック
## git clean -dfX ./ # 実行: 本当に消えるので実行するときは要注意

Introduction_to_Lattice_Boltzmann_Equation

目標

適当なタイミングで現代数学探険隊のプログラミング編を作りたい. そのための勉強.

メモ

  • 非圧縮性粘性流体の数値計算法.
  • 非圧縮粘性流体の特徴: 圧力を求める発展方程式がなく, 連続の式をみたす圧力をどう求めるか?
  • 歴史的には各時刻で圧力のポアソン方程式を解く.
    • これに時間がかかる.
  • 格子ボルツマン法 (Lattice Boltzmann Method, LBM)
    • 流体を微視的立場から捉える分子運動論を基礎にする
    • 格子気体モデル: 流体を速度をもつ有限個の仮想粒子の集合で近似する
    • 各粒子の衝突と並進を粒子の速度分布関数で逐次計算する
    • 速度分布関数のモーメントから巨視的流れ場を求める
    • 複雑な流れ場に対してもアルゴリズムが簡単
    • 質量・運動量の保存性がいい
    • 並列計算に適している

TODO

sample1.py

fem (FEM, 有限要素法)

参考メモ

rust

参考ページ

ode_dynamical_system001.rs

このページを参考にした. 可視化は python で gnuplot を生成してその gnuplot から png 生成し, ffmpeg で mp4 化. 後半にいくにつれて 10MB 程度ある csv を読み込んでは処理することになり, 処理が重くなるので軽量化したい.

ode_dynamical_system003.rs

ode_dynamical_system001.rs 同様. これは「回転」.

ode_dynamical_system004.rs

ode_dynamical_system001.rs 同様. これは「鞍点」. 可視化で新しいファイルを作らず, ode_dynamical_system003_visualize.py を転用している.

ode_dynamical_system005.rs

これは「不安定節」. 可視化は Python に切り替えた. いまの出力データ法からすると Gnuplot より Matplotlib を使う方がいい模様.

ode_dynamical_system006.rs

これは「安定節」.

ode_dynamical_system007.rs, ode_dynamical_system008.rs

うまく参考サイトのような図が描けなかった. 何が悪いのだろう?

ode_van_der_pol.rs

可視化は Gnuplot でやろうと思ったが, 方程式の解とベクトル場を同時に表示させられず断念した. Gnuplot の方が速そうなので, Gnuplot で何とかできるようになりたい.

ode_fitzhugh_nagumo.rs

  • TODO: 動画にパラメータの情報を入れて何をどう変えているかわかるようにする
  • YouTube へのリンク
  • math_memo.lyx に式と離散化メモあり
  • 予備考察として python 版を実行する

wave_eq

URL メモ

ファイルごとの方程式の説明

1dim_wave_eq_only_u.rs

このページを参考にして, 素直に波動方程式を解いている.

初期条件は $u = 0, u_t = 0$ で, 原点で波を $(1/2) \sin (2 \pi f t)$ で駆動していて, 離散化については math_memo.lyx を見ること.

調整中 1dim_wave_eq_only_u_with_src.rs

1dim_wave_eq_only_u.rs では u を直接動かしているが, こちらは外力 f で駆動しようとしている. 振幅が異様に小さく, なぜうまくいかないかまだわかっていない.

1dim_wave_eq_pml_both_side.rs

1dim_wave_eq_only_u.rs と同じ条件で波を生成している. 両端に PML をつけていて反射が起きない.

1dim_wave_eq_pml_both_side_compare.rs

1dim_wave_eq_pml_both_side.rs に対して PML なしの波も追加した. PML の有無での挙動を比較している.

1dim_wave_eq_pml_processing.rs

Qiita: Processingでシミュレーション~境界付近で波を消すにある processing のコード, processing/sketch_1dim_wave_eq_pml.pde を直接移植したコード. 左端で駆動していて領域中央の cnf.nx / 2 から吸収壁が始まる.

元ネタはおそらく次の英語のノート.

PML は変則的な 1 階化を使う. これは math_memo.lyx にまとめた. 勘違いしないようにドキュメントをよく読んで注意すること.

1dim_wave_eq_pml_right_wall.rs

1dim_wave_eq_pml_processing.rs と原則同じ. こちらは右端の 5 層の小領域で一気に吸収させている.

1dim_wave_eq_u_ut.rs

初期条件や波の駆動は 1dim_wave_eq_only_u.rs と同じ設定にしている. 2 階の常微分方程式でよくやるように, $u, v = u_t$ として, 次の連立の偏微分方程式を解く形にしている.

\begin{align} u_t = v, \quad v_t = u_{xx}. \end{align}

1dim_wave_eq_u_ut_with_src.rs

ふつうの 1 階化方程式で波源を中央につけた.

1dim_wave_eq_u_ut_with_src_pml.rs

波源と PML を両方つけた. 比較用に PML なしの近似解も計算している.

math_memo.lyx にある離散化の注意を守りつつ, 内部領域と PML 領域で解く方程式を変えている. 定式化については math_memo.lyx を見ること.

1dim_wave_eq_u_ut_with_src_pml_nonzero_bdval.rs

境界値を 0 からずらしてうまくいくかどうか検証した.

2dim_wave_eq_only_u.rs

中心を強制振動させ, それで波を起こしている. これだけだとわかりづらいが, PML つきの動画と比較すると反射が起こっていることがわかる.

2dim_wave_eq_pml_both_side.rs

上の動画と同じく中心を強制振動させ, それで波を起こしている. こちらは吸収壁をつけた. 上の動画と比較すると境界で反射が起きていないことがわかる.

PML の定式化については math_memo.lyx 参照.

TODO

次のリンク先の Julia コードを実装する.

Rust でもベクトル計算したいのだが, できる?

TODO リスト・参考文献

Python または Rust でコードを書いて貯めたい.

YouTube 投稿用メモ

追加すべきリスト

タイトル用サンプル

  • Rust 有限体積法 1 次元の線型移流方程式 アニメーションサンプル 数値流体解析の基礎 Visual C++とgnuplotによる圧縮性・非圧縮性流体解析

動画コメント用サンプル

コードは GitHub に置いてあるので興味があればどうぞ.

対応するのは次のコードです. ファイル・ディレクトリ整理で指定の場所にない場合は上の mathcodes から適当に探してください.

これ以外に, 数学・物理・プログラミングに関する無料の通信講座を運営しています. もしあなたがご興味あるなら, ぜひ次の通信講座ページ一覧を眺めてみてください.

他にも YouTube で数値実験の動画を投稿しています. リストにしているのでこちらもぜひどうぞ.

GitHub でもいくつかコードを公開しています.

ぜひチャンネル登録もしてください。

About

my codes for math or phys, and so on

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published