さていよいよ移流方程式を解くことを考えよう。以下では簡単のため、1次元の移流方程式
$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}=0$$
を解くことを考える。これは速度$u$で何らかの物理量$T$を形を変えずに運んでいることに相当する。
これはわざわざ数値計算などで解くまでもなく、右（$x$についてプラス側に動くものは$T(x,t)=f(x-ut)$という一般解を持つことがわかる。
さて問題は、これをどのようにして数値的に解くことができるか、ということである。これを差分形式で書いてみよう。
$$\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)$$
と考えるかもしれない。
まずこれを数値的に解くことにして、以下の差分形式でサンプルコードを作成してみよう
練習問題9-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");
}
}
}これは、実は数値的に発散してしまうことを確かめるためのコードになっている。これでは数値的に発散してしまうことが確かめられたら、次に差分の方法をやりなおしてみる。まず
$$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)$$
としてみよう。計算を100ステップごとに結果を書きだすように変更して、さらに時間時間を10から500に変更してたサンプルコードを完成してみよ：
練習問題9-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],cfl;
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=1; i<imax; i++)
T[i]=ここに記述する;
if (t%100==0){
for (i=1; i<imax; i++)
printf("%2.2e,",T[i]);
printf("\n");
}
}
}