さていよいよ移流方程式を解くことを考えよう。以下では簡単のため、1次元の移流方程式
$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}=0$$
を解くことを考える。これは速度$u$で何らかの物理量$T$を形を変えずに運んでいることに相当する。
これはわざわざ数値計算などで解くまでもなく、右（$x$についてプラス側)に動くものは$T(x,t)=f(x-ut)$という一般解を持つことがわかる。
さて問題は、これをどのようにして数値的に解くことができるか、ということである。これを差分形式で書いてみよう。物理量$T$について、時間を$n$、x座標上の位置を$i$という添え字で表すことにする。つまり時間$n$の時の場所$i$の物理量を$T^n_i$と書くことにすると、
$$\frac{T_i^{n+1}-T_i^{n}}{\Delta t}=-u\frac{T_{i+1}^n-T_i^n}{\Delta x}$$
$$T_i^{n+1}=T_i^{n}-u\frac{\Delta t}{\Delta x}(T_{i+1}^n-T_i^n)$$
と考えるかもしれない。
まずこれを数値的に解くことにして、以下の差分形式でサンプルコードを作成してみよう。例によって、paiza.io, onlineGDB, programiz, codechef.com, ideone.com, paiza.io などを利用して実行してみよう：
練習問題7-1.
#include <stdio.h>
#include <math.h>
#define u 1.0
#define dt 0.1
#define xmax 200
#define imax 200
#define tmax 10
int main(void){
double dx=(double)xmax/imax;
double x;
double T[xmax+1],cfl;
int t;
cfl=u*dt/dx;
int i;
for (i=1; i<=imax; i++){
if (i>=10&&i<50)
T[i]=1.0;
else
T[i]=0.0;
}
for (i=1; i<imax; i++)
printf("%2.2e,",T[i]);
printf("\n");
for (t=1; t<=tmax; t++){
for (i=1; i<imax; i++)
T[i]=ここに記述する;
if (t%1==0){
for (i=1; i<imax; i++)
printf("%2.2e,",T[i]);
printf("\n");
}
}
}これは、実は数値的に発散してしまうことを確かめるためのコードになっている。なお、ここで出力された値をグラフにするには、gnuplotのようなグラフを書くための便利な手法を使ってもらうのも良いし、MathematicaやMATLABなどを利用する手もあるが、excelやgoogle スプレッドシートでも十分に機能する。excelで描画するには、計算結果をコピペしてエクセルシートに貼り付けてみてから、「データ」タブから「区切り位置」を選んで、カンマ区切りにするという作業が必要な場合もある。計算結果をメモ帳などに一旦コピー・ペーストして、〇〇.csvなどとして保存し、これをexcelで読みこむ方法もある（csvは、カンマ区切りの数値データという意味で使われる拡張子）。
さて、上の方法では数値的に発散してしまうことが確かめられたら、次に差分の方法をやりなおしてみよう。まず
$$T_i^{n+1}=T_i^{n}-u\frac{\Delta t}{\Delta x}(T_{i+1}^n-T_i^n)$$
の代わりに
$$T_i^{n+1}=T_i^{n}-u\frac{\Delta t}{\Delta x}(T_{i}^n-T_{i-1}^n)$$
としてみよう。例によって、paiza.io, onlineGDB, programiz, codechef.com, ideone.com, paiza.ioなどを利用して、計算を100ステップごとに結果を書きだすように変更して、さらに時間時間を10から500に変更してたサンプルコードを完成してみよ：
練習問題7-2.
#include <stdio.h>
#include <math.h>
#define u 1.0
#define dt 0.1
#define xmax 200
#define imax 200
#define tmax 500
int main(void){
double dx=(double)xmax/imax;
double x;
double T[xmax+1],Tcopy[xmax+1],cfl; //Tcopyはi-1段階のものを入れることにする
int t;
cfl=u*dt/dx;
int i;
for (i=0; i<=imax; i++){
if (i>=10&&i<50)
T[i]=1.0;
else
T[i]=0.0;
}
for (i=1; i<imax; i++)
printf("%2.2e,",T[i]);
printf("\n");
for (t=1; t<=tmax; t++){
for (i=0; i<imax; i++)
Tcopy[i]=T[i]; // 前の値をTcopyにコピー
for (i=1; i<imax; i++)
T[i]=Tcopy[i]-cfl*(Tcopy[i]-Tcopy[i-1]); // 前のステップの結果 Tcopy を用いて計算
if (t%100==0){
for (i=1; i<imax; i++)
printf("%2.2e,",T[i]);
printf("\n");
}
}
}この結果は、
となって、かなり安定した計算ができていることに気が付く。しかし本来は青色の角ばった形状であったものが、500ステップ後には丸まった形へと変化してしまっている。これは特に角の部分の値が大きく変化するところが、次第になまっていくことに原因があることは容易に想像できる。空間の離散化をもっと細かくすると、少しは改善される部分があろうが、解決しないことも想像できるだろう。
いずれにせよ、風上側の差分を使って解く方が安定していることがわかった。これを風上差分と呼ぶ。しかし速度$u$が逆方向を向くと、差分表記から変更しないと発散してしまう可能性があるのは、都合が悪い。特に複雑な現象で、右にも左にも移動することがあるケースは、うまく解けないことになる。
課題：差分の方法をやりなおして、$$T_i^{n+1}=T_i^{n}-u\frac{\Delta t}{2\Delta x}(T_{i+1}^n-T_{i-1})$$としてみるとどうだろうか？たとえば
のような感じの結果（前進差分よりはマシだが、振動しながら進むという結果）などが得られると思う。与える条件を適宜考えて、こうした状況を確認し把握したことを簡単にまとめよ。
以上のように、少し差分形式を変えるだけで、得られる解が大きく変わってしまうことがある。これは差分法だけの問題ではないが、収束計算を行わない場合には特に問題となる場合がある。なお移流項を含むこの問題は、多項式で近似するなどして、特に変化の大きい部分について取り扱う手法（CIP法など）を用いて解決する場合が多い。
課題は講義時間内に解けるだけ解いてもらうことを想定しています。そのためできたところまでで良いので、講義時間中に提出してくれて構いません。自分で凝ったことをやってみたい人のために、締切を週末（11/2、深夜）としておきます。以下で提出してください。なお、g.eccアカウントでのアクセスが必要です。
課題の提出先はこちら：https://forms.gle/UNNyV5SEbj2Dv9FC7