2025年5月21日水曜日

記事一覧

all.md

数値計算の備忘録

数値計算で忘れそうなことを不定期でメモします。

PCからのほうが見やすいですが、スマホの場合は横画面にすると見やすいです。

<数値計算>

  1. オイラー法
  2. ルンゲクッタ法
  3. 逆行列の計算
  4. LAPACKのDGESVで連立方程式を解く
  5. 一様乱数生成 (fortran)
  6. マクスウェル分布に従う乱数生成 (fortran)
  7. FORTRANの関数のメモ
  8. FORTRANの行列・配列関連のメモ
  9. Gnuplotでカーブフィッティング
  10. Minpackで関数フィッティング (fortran)
  11. Quadpack で数値積分 (fortran)
  12. gslライブラリを使う (fortran)
  13. FFTW3 でフーリエ変換 (fortran)
  14. odepackで常微分方程式を解く (fortran)
  15. polyroots-fortran で多項式を解く (fortran)

<プラズマ関連>

  1. マクスウェル分布
  2. 単位 eV で表したマクスウェル分布
  3. 単位 eV のEEDFを m/s に変換する

<原子分子過程>

  1. 導入
  2. 断面積の取得
  3. 反応速度係数の計算

<その他・雑記事>

  1. LatexのBeamerでパワポ
  2. 原子質量単位のメモ
  3. ビットとバイト
  4. ポアソン方程式
  5. pythonのsliderで時系列データを可視化
  6. pyplot-fortranで fortran から直接グラフ生成

プログラムは FORTRAN で、Ubuntu OS のパソコンでやっています。

このサイトは Google Blogger で、Markdown で書いたものをHTML変換して作っています。

何か連絡がありましたらよろしくお願いします。

polyroots-fortran で多項式を解く

doc.md 多項式のお手軽なソルバー。
こういった複素数関連のソルバを自力で実装すると複素数の丸め誤差が出てしまうのでかなりありがたいところです


setup

前準備として fpm をインストールしておく (ココを参照).

githubからソースコードを持ってきてビルド

git clone https://github.com/jacobwilliams/polyroots-fortran.git cd polyroots-fortran/ fpm build --profile release fpm install sudo cp ~/.local/lib/libpolyroots-fortran.a /usr/local/lib/ sudo cp ~/.local/include/polyroots_module.mod /usr/local/include/ ls /usr/local/lib ls /usr/local/include

計算

yahoo知恵袋に 超難問 自作問題 があるのでこれを解いてみる。

32x^5 + 16x^4 -32x^3 -12x^2 + 6x + 1 = 0

program test use polyroots_module implicit none integer, parameter :: deg = 5 real(8) :: p(deg+1) = [32,16,-32,-12,6,1], zr(deg), zi(deg) integer :: i, istat call polyroots( deg, p, zr, zi, istat ) do i = 1, deg if(istat==0) write(*,'(2ES16.5)') zr(i), zi(i) end do end program test
  • deg で次数を指定
  • zrはリアルパート、ziはイマジナリーパート
  • p に係数を入力
  • call polyroots で計算
  • istat = 0 で成功

status output:

  • 0 : success
  • -1 : leading coefficient is zero
  • -2 : no roots found
  • >0 : the number of zeros found
gfortran test.f90 -L/usr/local/lib -lpolyroots-fortran -llapack -lblas -I/usr/local/include

lapackとblasもリンクしておく (toml参照)

$ ./a.out 8.41254E-01 0.00000E+00 -9.59493E-01 0.00000E+00 -6.54861E-01 0.00000E+00 4.15415E-01 0.00000E+00 -1.42315E-01 0.00000E+00

どうやら全て実数解になるらしい。 解は知恵袋の回答によると

x 〜 -0.959493 x 〜 -0.654861 x 〜 -0.142315 x 〜 0.415415 x 〜 0.841254

おそらくいい感じに計算できた。

2025年4月6日日曜日

OdePack で常微分方程式を解く

doc.md

インストール

準備として fpm をあらかじめインストールしておく (参考)。

git clone https://github.com/jacobwilliams/odepack.git cd odepack/src/ make clean make gfortran make run make test

自分の環境だと fpm build ではコンパイルできなかったので、Makefileを使った。
(external *** のように外部宣言している部分をコメントアウトすればうまくいくかも)
コンパイルすると M_odepack.o と m_odepack.mod が生成されるので、それぞれ /usr/local/lib、/usr/local/include などに移動しておく(任意).

計算準備

次の微分方程式を解く: $$ y'' -y' -2y = 18xe^{2x} $$

初期条件は $$ y(0) = \pi+\sqrt{2} = 4.5558062159628880 $$ $$ y'(0) = -2 + 2\pi - \sqrt{2} = 2.8689717448064913$$

理論解は以下の通り: $$ y' = (6x-2)e^{2x} + 2(3x^2-2x+\pi)e^{2x} - \sqrt{2}e^{-x} $$ $$ y = (3x^2-2x+\pi)e^{2x} + \sqrt{2}e^{-x} $$

ソースコード・コンパイル

lsoda を使って計算する。

program test use m_odepack implicit none ! odepack parameters integer, parameter :: neq = 2 real(8) :: atol(neq)=1d-9, rtol=1d-9 real(8) :: y(neq), t, tout, rwork(200) integer :: iopt=0 integer :: iout integer :: istate=1, itask=1 integer :: itol=2 integer :: jt = 2 integer :: liw=100, lrw=200, iwork(100) ! function setting integer, parameter :: N = 50 real(8), parameter :: tmin=0.d0, tmax=3.d0 real(8), parameter :: dt=(tmax-tmin)/dble(N-1) real(8) :: time(N) integer :: i t = tmin ! initial time y(1) = 4.5558062159628880d0 ! y(0) y(2) = 2.8689717448064913d0 ! y'(0) do i = 1, N if(i==1) time(i) = tmin if(i/=1) time(i) = time(i-1) + dt end do do i = 1, N tout = time(i) call dlsoda(fex,[neq],y,t,tout,itol,[rtol],atol,itask,istate,iopt,rwork,lrw,iwork,liw,jdum,jt) write (55,'(7ES12.4)') t, y(1), ftrue(t), y(1)-ftrue(t), y(2), dftrue(t), y(2)-dftrue(t) if ( istate<0 ) stop "istate < 0" enddo contains ! ================================================== ! subroutine fex( Neq, T, Y, Ydot ) integer, intent(in) :: neq real(8), intent(in) :: t, y(neq) real(8), intent(out) :: ydot(neq) Ydot(1) = Y(2) Ydot(2) = Y(2) + 2.d0*Y(1) + 18.d0*t*exp(2.d0*T) end subroutine fex ! ================================================== ! subroutine jdum() implicit none end subroutine jdum ! ================================================== ! real(8) function ftrue(x) real(8), intent(in) :: x real(8) :: k1 = sqrt(2.d0), k2 = acos(-1.d0) ftrue = (3.d0*x**2 -2.d0*x + k2)*exp(2.d0*x) + k1*exp(-x) end function ftrue ! ================================================== ! real(8) function dftrue(x) real(8), intent(in) :: x real(8) :: k1, k2 k1 = sqrt(2.d0) k2 = acos(-1.d0) dftrue = (6.d0*x - 2.d0)*exp(2.d0*x) + 2.d0*(3.d0*x**2 - 2.d0*x + k2)*exp(2.d0*x) - k1*exp(-x) end function dftrue ! ================================================== ! end program test
  • neq は微分方程式の階数 (変数の数)

  • atol は絶対誤差 (abs)、rtolは相対誤差 (rel)

  • その他は odepack の手続き変数

  • 解く時間は自分で設定した. 最初に t = tmin として、ループでは t -->> tout になるように時間が更新される。

  • $y(1)$ を 求める $y$、$y(2)$ を $y'$ とすると、$ydot(1)=y' = y(2)$、$ydot(2)=y(2)+2y(1)+18xe^{2x}$ と書けるので、これを解く。dotは微分を意味。

  • $x<0$ も解けるが、初期値が必要. 今は理論解がわかっているのですぐに試せる.

コンパイル

gfortran /usr/local/lib/M_odepack.o test.f90 #-I/usr/local/include/

実行結果

$ head fort.55 0.0000E+00 4.5558E+00 4.5558E+00 0.0000E+00 2.8690E+00 2.8690E+00 0.0000E+00 6.1224E-02 4.7554E+00 4.7554E+00 1.8396E-09 3.6747E+00 3.6747E+00 5.7735E-09 1.2245E-01 5.0092E+00 5.0092E+00 2.3605E-09 4.6483E+00 4.6483E+00 6.6694E-09 1.8367E-01 5.3288E+00 5.3288E+00 3.5579E-09 5.8302E+00 5.8302E+00 9.3923E-09 2.4490E-01 5.7283E+00 5.7283E+00 6.7996E-09 7.2697E+00 7.2697E+00 1.7355E-08 3.0612E-01 6.2254E+00 6.2254E+00 6.3543E-09 9.0258E+00 9.0258E+00 1.5290E-08 3.6735E-01 6.8415E+00 6.8415E+00 7.4869E-09 1.1170E+01 1.1170E+01 1.7392E-08 4.2857E-01 7.6028E+00 7.6028E+00 9.2457E-09 1.3788E+01 1.3788E+01 2.1080E-08 4.8980E-01 8.5416E+00 8.5416E+00 1.1599E-08 1.6984E+01 1.6984E+01 2.6115E-08 5.5102E-01 9.6968E+00 9.6968E+00 1.6157E-08 2.0880E+01 2.0880E+01 3.6313E-08 時刻 t 計算した y 理論解 yt 誤差 y-yt 計算したy' 理論解 y't 誤差 y'-y't

誤差は自分で設定した tol = 1d-9 のオーダーになっている。
例えば変数を減らす (ex: y' = ky) ときは neq = 1、増やすときは neq = 3 のようにする

参考

chatgptさんによると、

lsoda という解法について教えて! ChatGPT: いい質問です!LSODA は 常微分方程式 (ODE) を数値的に解くための非常に有名で賢いアルゴリズムのひとつです。 Fortran で最初に実装されましたが、今では Python(SciPy)、MATLAB、Julia など色々な言語で使えるようになっています。 🔧 LSODA の特徴 LSODA は 自動的に「剛性のある方程式かどうか」を判断して、適切な手法を切り替えてくれるハイブリッドなソルバ です。 剛性でない場合(stiffでない): → Adams method(マルチステップ法)を使用 剛性がある場合(stiff): → BDF法(後退差分公式) を使用 これによって、ユーザーが「剛性かどうか」を気にせずに使える、便利なソルバになっています。 📌 LSODA の使いどころ 方程式が剛性かどうか分からないとき 硬い部分とそうでない部分が混在しているような現象(例:化学反応、プラズマ、流体など) 長時間の安定したシミュレーションが必要なとき 🔍 Python(SciPy)での使い方例 python コピーする 編集する from scipy.integrate import solve_ivp import numpy as np import matplotlib.pyplot as plt # 例:dy/dt = -y, 初期値 y(0) = 1 def dydt(t, y): return -y sol = solve_ivp(dydt, [0, 10], [1], method='LSODA') plt.plot(sol.t, sol.y[0]) plt.xlabel("t") plt.ylabel("y") plt.title("dy/dt = -y を LSODA で解いた結果") plt.grid() plt.show() 🧠 補足:剛性(Stiffness)って? 剛性のある方程式とは、急激に変化する成分と緩やかに変化する成分が混在していて、普通の数値解法だとすごく小さなステップが必要になってしまうような ODE のこと。 LSODA のような「剛性対応ソルバ」がないと、解くのに時間がかかりすぎたり、不安定になったりします。 もし Fortran で使いたい場合や、実際に解きたい方程式があるなら、それに合わせてコード例を出すこともできるよ!どうする?

だそうです。。

2025年4月5日土曜日

pyplot-fortranで fortran から直接グラフ生成

doc.md

fortran のスクリプトから直接画像を生成できる。

インストール

fpm をあらかじめインストールしておく (参考).

git clone https://github.com/jacobwilliams/pyplot-fortran.git cd pyplot-fortran/ fpm build fpm test fpm install sudo mv ~/.local/lib/libpyplot-fortran.a /usr/local/lib sudo mv ~/.local/include/pyplot_module.mod /usr/local/include

ソースコード・コンパイル

program test use pyplot_module, only : pyplot, wp => pyplot_wp implicit none integer, parameter :: N = 100 real(8), parameter :: pi = 4.d0*atan(1.d0) real(8), parameter :: dx = 2.d0*pi / dble(N-1) real(8) :: x(N)=0.d0, y(N)=0.d0 integer :: i type(pyplot) :: plt integer :: istat ! set function do i = 2, N x(i) = x(i-1) + dx end do y = sin(x) ! plot graph call plt%initialize(grid=.true., xlabel='x', ylabel='y', figsize=[8,5], title='sinwave') call plt%add_plot(x,y, label='sin', linestyle='b-o', markersize=1, linewidth=2, istat=istat) call plt%savefig('graph.png', istat=istat) end program test
gfortran test.f90 -I/usr/local/include -L/usr/local/library/ -lpyplot-fortran ./a.out

結果

結構便利かも。
他にも build したときの test をみると pyplot でプロットできるグラフなら何でもプロットできそう。
出力過程は tmp.py に関数の数値データを全て入力した後に実行しているから、重いデータのプロットには時間がかかるかも

python の slider で時系列データを可視化

doc.md

インストール

pip install numpy pip install matplotlib

ボールの放物線運動

fortranで計算しておいた。

program ball real(8) :: t=0.d0, x=0.d0, y=0.d0, g=9.8d0 real(8) :: vx=3.d0, vy=50.d0, dt=0.2d0 do i = 1, 51 write(100,'(3f16.4)') t, x, y x = x + vx*dt y = y + vy*dt vx = vx vy = vy - g*dt t = t + dt end do end program ball

プロット

import numpy as np import matplotlib.pyplot as plt from matplotlib.widgets import Slider data = np.loadtxt('fort.100') time = data[:, 0] x = data[:, 1] y = data[:, 2] fig, ax = plt.subplots(figsize=(10, 5)) plt.subplots_adjust(bottom=0.25) line, = ax.plot([], [], marker='o', markersize=5, label='') ax.set_xlim(min(x), max(x)) ax.set_ylim(min(y), max(y)) ax.grid() ax_step = plt.axes([0.25, 0.1, 0.65, 0.03]) step_slider = Slider(ax_step, 'time', 0, len(data) - 1, valinit=0, valstep=1) def update(val): step = int(step_slider.val) line.set_data(x[:step + 1].tolist(), y[:step + 1].tolist()) ax.set_title(f'time = {time[step]:.2f} s') fig.canvas.draw_idle() step_slider.on_changed(update) plt.show()

結果

スライダーをずらすとその時刻までの軌道をプロットしてくれる。
時系列データのプロットに便利かも。

FFTW3 でフーリエ変換

doc.md

インストール

本家HP (リンク)に移動して http: fftw-3.3.10.tar.gz をダウンロード、適当な場所に移動しておく。

tar -zxvf http: fftw-3.3.10.tar.gz cd fftw-3.3.10 ./configure --prefix=/usr/local/ make sudo make install

他にも並列化オプションなどが必要であれば設定しておく。

ソースコード・コンパイル

10Hz, 20Hz, 30Hzの sin 信号を足し合わせた以下の信号をFFTする:
$$f(t) = \sin(2\pi 10t) + 0.7\sin(2\pi 20t) + 0.3\sin(2\pi 30t)$$ もちろん、FFTスペクトルにはf=10, 20, 30 に y= 1, 0.7, 0.3 のピークが立つ.

program test implicit none include "fftw3.f" integer, parameter :: N = 100 real(8), parameter :: dt=1.d0/dble(N-1), fs=1.d0/dt real(8) :: t(N), freq(N), f(N) complex(8) :: cf(N) ! f(t) を複素数配列に変換したもの complex(8) :: f1(N) ! FFT[ f(t) ] complex(8) :: f2(N) ! INVFFT[ FFT[ f(t) ] ] integer(8) :: plan real(8), parameter :: K=2.d0*4.d0*atan(1.d0) integer :: i ! set f(t) do i = 1, N t(i) = real(i-1) * dt freq(i) = dble(i-1) * fs / dble(N-1) f(i) = sin(K*10.d0*t(i)) + 0.7d0*sin(K*20.d0*t(i)) + 0.3d0*sin(K*30.d0*t(i)) cf(i) = cmplx(f(i), 0.0d0) ! real -->> complex end do ! FFT[ f(t) ] call dfftw_plan_dft_1d(plan, N, cf, f1, FFTW_FORWARD, FFTW_ESTIMATE) call dfftw_execute_dft(plan, cf, f1) call dfftw_destroy_plan(plan) call dfftw_cleanup() ! inversion FFT call dfftw_plan_dft_1d(plan, N, f1, f2, FFTW_BACKWARD, FFTW_ESTIMATE) call dfftw_execute_dft(plan, f1, f2) call dfftw_destroy_plan(plan) call dfftw_cleanup() ! amplitude setting f1 = f1/dble(N/2) f2 = f2/dble(N) do i = 1, N write(55,'(3ES16.2)') t(i), f(i), real(f2(i)) write(56,'(4ES16.2)') freq(i), real(f1(i)), aimag(f1(i)), abs(f1(i)) end do end program test
gfortran test.f90 -I/usr/local/include/ -L/usr/local/library -lfftw3 ./a.out

結果

問題なく計算できた。

エイリアシングも見えている。

その他

program test implicie none include "fftw3.f"

のように fftw3.f (変数定義のみされたファイル) を include する方式は使いにくい。
できれば use fftw3 のようにしたいところ。
この場合、fftw3.f の先頭行に module fftw3、末行に end module fftw3 と書いて、別途コンパイルしておいてもよいかも。

fortran で gsl を使う

doc.md

インストール

sudo apt install libgsl-dev

ソースコード・コンパイル

program bessel use iso_c_binding implicit none interface function gsl_sf_bessel_J0(x) bind (c, name="gsl_sf_bessel_J0") import :: c_double real(c_double) :: gsl_sf_bessel_J0 real(c_double), value, intent(in) :: x end function gsl_sf_bessel_J0 end interface print*, gsl_sf_bessel_J0(5.d0) end program bessel
gfortran test.f90 -lgsl #-lgslcblas ./a.out -0.17759677131433826 (ベッセル関数 J0(5.0)の値)

ベッセル関数は大文字の J と小文字の j で異なるから注意。
使える関数はマニュアルを見るとわかる。
ベクトル (配列) を扱う関数についてはテスト中 (うまくいってない)。