一次元移流方程式をルンゲ-クッタ法で時間高次精度化してC言語で解く

#陽的4段4次ルンゲクッタ法
ルンゲクッタ法は微分方程式の数値計算を行うための近似解を与える手法であり、今回は4次ルンゲクッタ法を採用した。その理由として、ルンゲクッタ法では4次の時までは計算コストに比例して計算精度が上がっていくが、4次以上では計算コストの増加と計算の正確さが見合わなくなっていくことから、一般に数値計算では4次以上のルンゲクッタ法は使用されないためである。

ここで、陽的4段4次ルンゲクッタ法の計算スキームは次式で与えられる。ただし空間差分は2次の中心差分を採用しており、$Δt$は時間刻み幅、下付きの$i$は空間方向、上付きの$n$は時間方向の格子点、$h$はステップ幅を表す。

x_{n+1}=x_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4) \tag{1.10}

ここで、$k_1,k_2,k_3,k_4$は次式で与えられる。

\begin{align}
&k_1=hf(x_n,t_n) \\
&k_2=hf(x_n+\frac{1}{2}k_1,t_n+\frac{1}{2}h) \\
&k_3=hf(x_n+\frac{1}{2}k_2,t_n+\frac{1}{2}h) \\
&k_4=hf(x_n+k_3,t_n+h) \\
\end{align}

#ルンゲクッタ法による一次元移流方程式の時間高次精度化
今回は、以下の以下の一次元移流方程式を解く。

\frac{df}{dt}+c\frac{df}{dx}=0 \tag{1.11}

一次元移流方程式(1.11)をルンゲクッタ法(1.10)を用いて風上差分で離散化すると、次式(1.12)となる。ただし、$Δt$は時間刻み幅、下付きの$i$は空間、上付きの$n$は時間の格子点を表し、風上側は常に$f_{i-1}$とする。

f^{n+1}_{i}=f^n_i-\frac{1}{6}\frac{\Delta t}{\Delta x}c(k_1+2k_2+2k_3+k_4) \tag{1.12}

ただし、$k_1,k_2,k_3,k_4$は次式で与えられる。

\begin{align}
&k_1=f^{n-1}_{i}-f^{n-1}_{i-1} \\
&k_2=f^{*}_{i}-f^{*}_{i-1} \\
&k_3=f^{**}_{i}-f^{**}_{i-1} \\
&k_4=f^{***}_{i}-f^{***}_{i-1} \\
\end{align}

また、上式の各アスタリスク$^*$つきの関数は次式で与えられる。

\begin{align}
&f^{*}_{i}=f^{n-1}_{i}+\frac{1}{2}k_{1}(i) \\
&f^{**}_{i}=f^{n-1}_{i}+\frac{1}{2}k_{2}(i) \\
&f^{***}_{i}=f^{n-1}_{i}+k_{3}(i) \\
&f^{*}_{i-1}=f^{n-1}_{i-1}+\frac{1}{2}k_{1}(i-1) \\
&f^{**}_{i-1}=f^{n-1}_{i-1}+\frac{1}{2}k_{2}(i-1) \\
&f^{***}_{i-1}=f^{n-1}_{i-1}+k_{3}(i-1) \\
\end{align}

式の導出は下記のページを参考にしました。
[https://www.ep.sci.hokudai.ac.jp/~gfdlab/comptech/resume/270_adveq3-1/2011_0303-takayasu.pdf]
(https://www.ep.sci.hokudai.ac.jp/~gfdlab/comptech/resume/270_adveq3-1/2011_0303-takayasu.pdf)

#解析条件
解析条件として、一次元移流方程式(1.11)の係数$c$は$1.0$とした。空間の刻み幅は $Δx=0.01$、格子点数 $i_{max}=100$ 固定として、繰り返しタイムステップ数 $n_{max}$は任意の値を設定し、時間の刻み幅 $Δt=1.0/n_{max}$とした。
また、初期波形としてFig.1および(1.21)式に示す振幅が$1.0$の正弦パルスを初期条件として与えた。

f(x) = \left\{
\begin{array}{ll}
\frac{1}{2} (cos(\pi+\frac{x-0.1}{0.1}\pi)+1) & (0.1 \lt x \lt 0.3) \\
0 & (x \leq 0.1, x \geq 0.3)
\end{array}
\right. \tag{1.21}

Fig.1 初期波形

#プログラム実装
プログラムはC言語で作成した。また、周期境界条件となるように実装を行った。境界については、始点と終点の外側にそれぞれ境界外側の格子点を5点取っている。これは、より高次の空間差分のため多めに外側の格子点を取っているためで、実際は今回の差分式では1点で良い。詳細については、下記のソースコードを参照願いたい。

/*********************************
     1D Advection Equation
  4th order accuracy Runge-Kutta
  GPL Copyright       2017/12/15
*********************************/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define Xnum 100
#define Nnum 400
#define Xmin 0.0
#define Xmax 1.0
#define dx ((Xmax-Xmin)/Xnum)
#define dt (1.0/Nnum)
#define NumBytes ((Xnum+11)*sizeof(double))

//Initial Condition
void Init(double *f, double *x, double *fold, double *init){
  FILE *fileout;
  int i,j,k;
  double xx;

  for(i=0;i<=Xnum+10;i++){
    x[i] = (double)(i-5)/(double)Xnum;
  }

  for(j=0;j<=Xnum+10;j++){
    if(x[j]>0.1 && x[j]<0.3){
      xx = (x[j]-0.1)/0.2;
      f[j] = 0.5*(cos(M_PI+2*M_PI*xx)+1);
    }
    else{
      f[j] = 0;
    }
  }

  for(j=0;j<=Xnum+10;j++){
    fold[j] = f[j];
  }

  fileout = fopen("output00000.dat","w");
  for(k=0;k<=Xnum+10;k++){
    fprintf(fileout,"%lf %lf\n",x[k],f[k]);
  }

}

void Diff(double *f, double *x, double *fnew, double *fold, double *init){
  FILE *fileout;
  char filename[99];
  int n,i,j,k;
  double dfdx,k1,k2,k3,k4,kn;

  for(n=1;n<=Nnum;n++){

    for(i=5+1;i<=Xnum+10;i++){
      //Rnnge-Kutta
      k1 = fold[i]-fold[i-1];
      k2 = (fold[i]+0.5*k1*x[i])-(fold[i-1]+0.5*k1*x[i-1]);
      k3 = (fold[i]+0.5*k2*x[i])-(fold[i-1]+0.5*k2*x[i-1]);
      k4 = (fold[i]+k3*x[i])-(fold[i-1]+k3*x[i-1]);
      kn = -1.0*(dt/dx)*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;

      fnew[i] = f[i] + kn;
    }

  //Save the old timestep
  for(j=0;j<=Xnum+10;j++){
    fold[j] = f[j];
  }

  //Update considering boundary condition
  for(j=5+1;j<=Xnum+5;j++){
    f[j] = fnew[j];
  }
  for(j=0;j<=5;j++){
    f[5-j] = fnew[Xnum+5-j];
  }
  for(j=1;j<=5;j++){
    f[Xnum+5+j] = fnew[5+j];
  }

  //Output
  sprintf(filename, "output%05d.dat",n);
  fileout = fopen(filename,"w");
  for(k=0;k<=Xnum+10;k++){
    fprintf(fileout,"%lf %lf\n",x[k],f[k]);
  }

  }
  printf("CFL number=%lf\n",dt/dx);
}


int main(void){
  double *f,*x,*fold,*fnew,*init;
  f    = (double *)malloc(NumBytes);
  fnew = (double *)malloc(NumBytes);
  fold = (double *)malloc(NumBytes);
  x    = (double *)malloc(NumBytes);
  init = (double *)malloc(NumBytes);

  Init(f,x,fold,init);
  Diff(f,x,fnew,fold,init);

  return 0;
}

#解析結果
はじめに、$Δx=0.01$、$Δt=0.005$、クーラン数 $C=0.5$とした時のオイラー法の計算結果をFig.2、4次ルンゲクッタ法の計算結果をFig.3にそれぞれ示す。なお、オイラー法の計算結果は比較用として示したが、計算式やソースコードなどについては割愛する。なお、クーラン数$C$の定義は次式(1.31)である。

クーラン数\  C=c\frac{\Delta t}{\Delta x} \tag{1.31}

Fig.2　クーラン数C=0.5　オイラー法  Fig.3　クーラン数C=0.5　4次ルンゲクッタ法

Fig.3より、クーラン数$C=0.5$の時の4次ルンゲクッタ法は解が発散してしまっていることが分かる。これは、時間高精度化した4次ルンゲクッタ法では、CFL条件を満たすため$C<0.5$でなければならず、今回の$C=0.5$ではこれを満たさないためである。従って、得られた結果は数値解析結果としては正しいと考えられる。


次にCFL条件を満たすように$Δx=0.01$、$Δt=0.0025$、クーラン数$C=0.25$とした時のオイラー法の計算結果をFig.4、4次ルンゲクッタ法の計算結果Fig.5 にそれぞれ示す。

Fig.4　クーラン数C=0.25　オイラー法  Fig.5　クーラン数C=0.25　4次ルンゲクッタ法

Fig.5より、クーラン数$C=0.25$の時の4次ルンゲクッタ法は発散せずに正しく計算できていることが確認できる。また、Fig.4のオイラー法の結果とFig.5を比較すると、初期波形の移流に伴う振幅の低下がFig.4のオイラー法よりも明らかに抑えられており、時間高精度化の効果が確認できる。従って、得られた結果は数値解析結果として正しいと考えられる。