はじめに

現在，私は研究で使用するため，熱のシミュレーションをC++を使って開発しています．しかし，C++は高速であるが，それ以外の利点は少ないと感じています．コードは冗長的で，デバッグも煩雑であるため，開発とは別のところで時間が奪われてしまうことが時々あります．またC++は描画ツールが乏しくカバーするためにPythonを使用しており，開発言語が2つあり，手間がかかると感じています．そんな時，先輩がJuliaは速く，書きやすいとJuliaの良さを熱弁していたため，どれくらいの速度がでるのかC++とJuliaで比較してみようと思います．

実行環境

OS: Ubuntu 24.04 (WSL2)
CPU: 16 × AMD Ryzen 7 5700X 8-Core Processor
RAM: 15 GiB

Juliaについて

Juliaについて簡単に説明します．
Juliaは，インタプリタ言語とコンパイル言語の特長を兼ね備えた，高性能なプログラミング言語です．動的型付けをサポートしつつ，コンパイル言語に匹敵する高速な実行速度を実現しています．この高速性を支える仕組みとして，JuliaはJust-In-Time (JIT) コンパイルを採用しています．JITコンパイルでは，コードが実行時にコンパイルされ，最適化されます．その結果，初回の実行には若干の時間がかかるものの，コンパイル後のキャッシュが利用されるため，2回目以降の実行は非常に高速です．この最適化により，JuliaはCやFortranといった従来のコンパイル言語に匹敵する処理速度を持つとされています．

コード

以下のレポジトリにあります．

比較の仕方

比較のために解く方程式は次の2次元の熱方程式とします．

$$
\frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
$$

$u(x, y, t)$は温度分布を表し、$\alpha$は熱伝導率を示します．

以下に偏微分方程式を離散化する方法について説明します．

離散化の方法

時間の離散化:
オイラー法 (Euler Method) を使用する．これは，時間ステップ $\Delta t$ ごとに温度分布を更新する方法です．

$$
u^{n+1} = u^n + \Delta t \cdot f(u^n)
$$

$u^n$ は時間ステップ $n$ における温度分布， $f(u^n)$は熱方程式の右辺です．

空間の離散化:
中心差分法 (Central Difference Method) を使用して空間の2次微分 $\frac{\partial^2 u}{\partial x^2}$ と $\frac{\partial^2 u}{\partial y^2}$を離散化します．

$$
\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \
, \
\frac{\partial^2 u}{\partial y^2} \approx \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{(\Delta y)^2}
$$

つまり，熱方程式の右辺である$f(u^n)$は次のようになります．

$$
f(u^n) = \alpha \left( \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{(\Delta y)^2} \right)
$$

よって，解くべき方程式は次のようになります．
$$
u_{i,j}^{n+1} = u_{i,j}^n + \Delta t \cdot \alpha \left( \frac{u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n}{(\Delta x)^2} + \frac{u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n}{(\Delta y)^2} \right)
$$

シミュレーションの設定

初期条件の設定
初期条件の温度は中心からの距離に応じて温度が指数関数的に減少するような分布を作成しています．
$r_0$は半径（75），$(x_c, y_c)$は空間の中心の座標，$(x,y)$は任意のセルの中心点を示します．

$$
u(x, y, 0) = 10 \cdot \exp\left(-\frac{\sqrt{(x - x_c)^2 + (y - y_c)^2}}{r_0}\right)
$$

境界条件の設定
境界の温度は0としています

メッシュの設定
空間を150×150と定義し，150×150のメッシュに分割して計算しています．

時間
シミュレーション時間は100秒行い，時間ステップ $\Delta t$は0.1で計算します．

安定解析

この方程式が発散しないことを示すために，CFL条件を求めます．線形方程式であるため，フォン・ノイマン安定性解析を使用します．フォン・ノイマン安定性解析では，解の形式を次のように仮定します．
$$
u_{i,j}^n = \hat{u}^n e^{ik_x x_i} e^{ik_y y_j}
$$

ここで，$\hat{u}^n$は振幅，$k_x$および$k_y$は波数です．この仮定を離散化された方程式に代入し，増幅因子を求めます．増幅因子は，解が時間とともにどのように振る舞うかを示します．

安定性パラメータ$r$は，$x$方向だけでなく$y$方向にも依存するため，次のように表現されます．
$$
r_x = \frac{\alpha \Delta t}{(\Delta x)^2}, \quad r_y = \frac{\alpha \Delta t}{(\Delta y)^2}
$$

そして，増幅因子$G$は次のように修正されます．
$$
G = 1 - 2r_x \cdot \left( \cos(k_x \Delta x) - 1 \right) - 2r_y \cdot \left( \cos(k_y \Delta y) - 1 \right)
$$

ここで，$r_x = \frac{\alpha \Delta t}{(\Delta x)^2}$は$x$方向の安定性パラメータ，$r_y = \frac{\alpha \Delta t}{(\Delta y)^2}$は$y$方向の安定性パラメータです．

安定性のためには，増幅因子の絶対値が1以下である必要があります．したがって，次の条件が求められます．
$$
|G| \leq 1
$$

この条件を満たすために，次のCFL条件が得られます．
$$
r_x + r_y \leq \frac{1}{2}
$$

つまり，$x$方向と$y$方向の両方に対する安定性パラメータを考慮することで，次のCFL条件が得られます．
$$
\frac{\alpha \Delta t}{(\Delta x)^2} + \frac{\alpha \Delta t}{(\Delta y)^2} \leq \frac{1}{2}
$$

この式により，$x$と$y$方向の空間ステップ$\Delta x$および$\Delta y$の両方に対して安定性条件を満たすことが必要であることが示されます．

$\Delta t=0.1$，それ以外はすべて1と設定しているため，条件を満たしているといえます．

比較方法

計測の仕方

計測はシミュレータを実行した箇所を測定しています．つまり初期化やファイル保存部分は測定外とします．計測は3回測定し，その平均を比較します．ただし，JuliaはJITコンパイルによって実行時間が長くなるため初回実行時ではなく，2回目以降の実行時間を使用します．
またC++の最適化フラグの有無についても調査してみます．最適化フラグは-O3 -march=native -fltoを使用します．（@fujitanozomuさん，コメントありがとうございます．）
以下はC++とJuliaの測定箇所です．

C++の測定箇所

  auto start_time = std::chrono::high_resolution_clock::now();
  for (int t = 0; t < num_time_steps; ++t) {
    solve_heat_equation(temperature, time_step, alpha, deltaX, deltaY);
  //  結果の保存関数
  //  save_to_file(temperature, t);
  }
  auto end_time = std::chrono::high_resolution_clock::now();

  std::chrono::duration<double> execution_time = end_time - start_time;
  std::cout << "Total execution time: " << execution_time.count() << " seconds"
            << std::endl;

Juliaの測定箇所

	start_time = now()
	solution = solve(problem, Euler(), dt = time_step, saveat = time_step)
  end_time = now()  

  println("Total execution time: $(end_time - start_time)")

結果

回数	C++[s] 最適化なし	Julia[s]	C++[s] 最適化あり
1	0.606	0.205	0.033
2	0.680	0.577	0.033
3	0.635	0.176	0.033
平均	0.640	0.319	0.033

コンパイルで最適化なしの場合，C++の実行時間平均は0.640秒で，Juliaと比較して，約2.0倍低速とJuliaのほうが高速でした．しかし，コンパイルで最適化ありだと，C++の実行時間平均は 0.033秒 で，Juliaと比較して約 9.7倍高速 になりました．

結論

C++のコンパイルが最適化される前はJuliaが約2.0倍高速でしたが，コンパイルを最適化すると，Juliaよりも9.7倍高速になり，C++がめちゃくちゃ高速であるとわかりました．ただ，それらしく記述するだけでJuliaも十分高速でしたので，使う機会があれば，使ってみたいと思います．これからもC++からは逃れそうにないです．今後もC++で頑張っていきます．
もしさらなる高速化や改善案があれば、ぜひ教えてください。

追記

記事投稿後，プログラミング言語の速度比較サイトってよくよく考えたら普通にありそうと思い，探してみたらあったので，備忘録のために残しておきます．

マンデルブロ集合を生成するコードを使用した速度比較
https://benchmarksgame-team.pages.debian.net/benchmarksgame/performance/mandelbrot.html
適応的台形法による数値積分アルゴリズムを用いた速度比較
https://github.com/astrojhgu/adaptrapezoid_benchmark