カテゴリー: Uncategorized

さていよいよ移流方程式を解くことを考えよう。以下では簡単のため、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");
    }
  }  
} 

れを差分形式で書いてみよう。

$$\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)$$

と考えるかもしれない。それで以下のようなコードでこれを検討してみよう。

#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]=T[i]-cfl*(T[i+1]-T[i]);
    if (t%1==0){
      for (i=1; i<imax; i++)
        printf("%2.2e,",T[i]);
      printf("\n");
    }
  }  
}

もともと、

のような形を初期条件としている。これが一定速度で右に動いて行ってもらいたいと考えているわけだが、実際には

となってしまい、計算が安定してくれない。

３次方程式 $2x^3+10x^2-2x=10$の解をニュートン法で求めてみる。$(x+1)(x-1)(2x+10)=0$とできることから解は3つあるが、ここでは正のものについて、初期値10からスタートしてみることにする。導関数 $6x^2+20x-2$ を利用すると、以下のようなコードを書くことができる。

#include <stdio.h>
#include <math.h>
#define eps 1.0e-15
#define xinit 10.0
double f(double x){
  return 2*x*x*x+10.0*x*x-2.0*x-10.0;
}
double fp(double x){
  return 6.0*x*x+20.0*x-2.0;
}
int main(void){
  double x0, x1, f0, fp0;
  x0=xinit;
  do {
    f0=f(x0);fp0=fp(x0);
    x1=x0-f0/fp0;
    x0=x1;
    printf(" %e %e \n", x0, f0);        
  }while(fabs(f0)>eps);
  printf("answer: %e \n", x0);
}

解の収束が悪かったり、似たような値の解が存在する場合などは、少し工夫が必要となる（減速ニュートン法やセカント法と呼ばれるものなど）。

（以下では例によって、開発環境が手元にない場合は、paiza.io, onlineGDB, programiz, codechef.com, ideone.com, paiza.ioなどを実行してみると良い）

連立常微分方程式も、同様の手法で解くことができる。たとえば生物の捕食―被食関係における個体数の変動を表現する数理モデルに、ロトカ・ヴォルテラの方程式がある。捕食者（たとえばヤマネコ）の個体数$x$が、被食者（たとえばウサギ）の個体数$y$とが、以下のような2元連立非線形常微分方程式の関係であらわされるとする考え方である：

$$\frac{dx}{dt}=ax – cxy$$

$$\frac{dy}{dt}=-by+cxy$$

こうした連立微分方程式であっても、同様にルンゲ・クッタ法を適用することで、現実的に解くことができる。やり方としては同じだが

$$\frac{dx}{dt}=f(t,x,y)$$

$$\frac{dy}{dt}=g(t,x,y)$$

と整理して、各々に1段ずつルンゲクッタ法を適用していけばよい。すなわち、

$$\begin{matrix}k_{x1}=hf(t_i,x_i,y_i) & k_{y1}=hg(t_i,x_i,y_i)\\k_{x2}=hf(t_i+\frac{h}{2},x_i+\frac{k_{x1}}{2},y_i+\frac{k_{y1}}{2})&k_{y2}=hg(t_i+\frac{h}{2},x_i+\frac{k_{x1}}{2},y_i+\frac{k_{y1}}{2}) \\k_{x3}=hf(t_i+\frac{h}{2},x_i+\frac{k_{x2}}{2},y_i+\frac{k_{y2}}{2})&k_{y3}=hg(t_i+\frac{h}{2},x_i+\frac{k_{x2}}{2},y_i+\frac{k_{y2}}{2})\\k_{x4}=hf(t_i+h,x_i+k_{x3},y_i+k_{y3}) &k_{y4}=hg(t_i+h,x_i+k_{x3},y_i+k_{y3})\end{matrix}$$

として順次計算し、

$$x_{i+1}=x_i+\frac{1}{6}(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})$$

$$y_{i+1}=y_i+\frac{1}{6}(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})$$

として次のステップの $x, y$を計算していけばよい。一見複雑に見えるが、やっていることは極めて単純で、これをプログラムにすると

#include <stdio.h>
#include <math.h>
#define N 1000
#define h 1.0
double f(double t,double x, double y);
double g(double t,double x, double y);
int main(void) {
  int i;
  double x=1000.0;
  double y=100.0;
  double t=0.0;
  double kx1,ky1,kx2,ky2,kx3,ky3,kx4,ky4;
  for (i = 0; i < N; i++) {
    kx1=（ここに記述する）; ky1=（ここに記述する）;
    kx2=（ここに記述する）; ky2=（ここに記述する）;
    kx3=（ここに記述する）; ky3=（ここに記述する）;
    kx4=（ここに記述する）; ky4=（ここに記述する）;
    x=（ここに記述する）
    y=（ここに記述する）
    t=t+h;
    printf("%4.4lf\t %4.4lf\t %e\n",x,y,t);
    }
}
double f(double t,double x, double y) {
  return (0.01*x-0.00001*x*y);
}
double g(double t,double x, double y) {
  return (-0.02*y+0.00001*x*y);
}

のように書くことができる。この結果が

1009.0456 99.0095 0
1018.1831 98.0376 1
1027.4132 97.0843 2
1036.7367 96.1491 3
1046.1546 95.2318 4
1055.6675 94.3323 5
1065.2765 93.4501 6
1074.9824 92.5852 7
1084.7859 91.7372 8
1094.6881 90.9059 9

のようになっていれば、計算はできている。1000回くらい繰り返すと、

のような餌と捕食者の変動の周期的な関係を見ることができる（モデルやパラメタの問題点は議論しない）。

これで常微分方程式に基づいた、さまざまな物理シミュレーションを実施するための基礎的な知識は身に付いた。なお、上記のコードは、ポインタの概念を使えば関数の読み出しの部分で断然簡素化することができるので、興味のある人はやってみよう。

今週の課題

課題は講義中にやってもらうことを想定しています、以下について回答してみてください。オリジナリティが高い回答ほど高評価で、おなじオリジナリティであれば提出時間が早いほど高評価。提出は10/10 深夜まで。提出が遅くてもオリジナリティが非常に高いものが最高評価、とします。提出はITCではなくて、https://forms.gle/5SAxcgPD1kkkKZ7w9 にお願いします。

ルンゲクッタの４次公式は、以下のようなものであった：

微分方程式$\frac{dy}{dx}=f(x,y)$ に対し、初期条件$y(x_0)=y_0$としたときの数値解は、刻み幅を$h$として、$x_n=x_0+nh$における$y$の値を$y_n$から$y_{n+1}=y_n+k$と書く場合、$k$は以下で近似できる
$$\begin{array}{ccc}k_1 & = & hf(x_n, y_n)\\k_2 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_1}{2}\right )\\k_3 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_2}{2}\right )\\k_4&=&hf(x_n+h,y_n+k_3)\\k&=&\frac{1}{6}(k_1+2k_2+2k_3+k_4) \end{array}$$

連立微分方程式
$$\left \{ \begin{array}{ccc}\frac{dy}{dx} &=&f(x,y,z),&初期条件が& y(x_0)=y_0\\\frac{dz}{dx} &=&g(x,y,z),&初期条件が& z(x_0)=z_0 \end{array} \right .$$
の数値解は、刻み幅を$h$、$x_n=x_0+nh$として以下で与えられる：
$$\begin{array}{ccc}k_1 & = & hf(x_n, y_n, z_n)\\l_1 & = & hg(x_n, y_n, z_n)\\k_2 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_1}{2},z_n+\frac{l_1}{2}\right )\\l_2 & = & hg \left ( x_n+\frac{h}{2},y_n+\frac{k_1}{2},z_n+\frac{l_1}{2}\right )\\k_3 & = & hf \left ( x_n+\frac{h}{2},y_n+\frac{k_2}{2},z_n+\frac{l_2}{2}\right )\\l_3 & = & hg \left ( x_n+\frac{h}{2},y_n+\frac{k_2}{2},z_n+\frac{l_2}{2}\right )\\k_4&=&hf(x_n+h,y_n+k_3, z_n+l_3)\\l_4&=&hg(x_n+h,y_n+k_3, z_n+l_3)\\k&=&\frac{1}{6}(k_1+2k_2+2k_3+k_4)\\l&=&\frac{1}{6}(l_1+2l_2+2l_3+l_4) \\y_{n+1}&=&y_n+k\\z_{n+1}&=&z_n+l\end{array} $$

上の考え方を使って、以下の問題を解いてみよう

課題1：$y^\prime =y-10x +1$を刻み幅0.1で$x=0$から$x=1.0$まで解け。ただし初期条件を$y(0)=1$とする。

課題２：2階微分方程式$\frac{d^2 y}{dx^2}-2x\frac{dy}{dx}+xy=0$を初期条件$y(0)=1, y^\prime(0)=1$で解け。

課題3: 野球ボールがどこまで飛ぶか考えよ。ただし空気抵抗が速度の2乗に比例する以下の形で書けるとする。
$$\begin{array}{ccc}a_x&=&-\frac{3}{4D}\frac{\rho_{air}}{\rho_{ball}} C_d V V_x\\a_y&=&-g-\frac{3}{4D}\frac{\rho_{air}}{\rho_{ball}}C_d V V_y \\V&=&\sqrt{V_x^2+V_y^2}\end{array}$$
ここで$a_x$を水平方向の加速度、$a_y$を垂直方向の加速度、$D$をボールの直径=7[cm]、$g$を重力加速度=9.8[m/s$^2$]、$\rho_{air}$を空気の密度=1.2[kg/m$^3$]、$\rho_{ball}$をボールの密度=650[kg/m$^3$]、$C_d$を抵抗係数=0.2とする。初速を50[m/s]、バットで打った場所が高さ1.0[m]である場合、45度で飛び出す時が最大距離ではないことを何らかの形で示せ。なおボールの回転による揚力や、風など空気の揺らぎは考えない。
なお、以下のコードを参考にしても良い。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 1000
#define h 0.01
double G=-9.8;
double k=-0.00396;
double f(double t,double x, double y);
double g(double t,double x, double y);
int main(void) {
  int i;
  double v=50.;
  double degree=38.;
  double radian=degree*3.141592/180.0;
  double x=v*cos(radian);
  double y=v*sin(radian);
  double t=0.0;
  double xx=0.,yy=1.0;
  double vx, vy,x1,x2,x3,x4,y1,y2,y3,y4;
  double kx1,ky1,kx2,ky2,kx3,ky3,kx4,ky4;
  for (i = 0; i < N; i++) {
    x1=ここに記述;y1=ここに記述;
    kx1=ここに記述;ここに記述;
    vx=ここに記述;ここに記述;
    x2=ここに記述;y2=ここに記述;
    kx2=ここに記述;ky2=ここに記述;
    vx=ここに記述;vy=ここに記述;
    x3=ここに記述;y3=ここに記述;                  
    kx3=ここに記述;ky3=ここに記述;
    vx=ここに記述;vy=ここに記述;
    x4=ここに記述;y4=ここに記述;
    kx4=ここに記述;ky4=ここに記述;
    xx=ここに記述;
    yy=ここに記述;    
    x=ここに記述;
    y=ここに記述;
    t=t+h;
    if (yy<0.0){
      printf("max=%4.1lfm\t%4.4lf\t %4.4lf\t %e\t %4.4e\t %4.4e\n",xx,x,y,t,xx,yy);
      exit(0);
    }
  }
}
double f(double t,double x, double y) {
  return (k*sqrt(x*x+y*y)*x);
}
double g(double t,double x, double y) {
  return (G+k*sqrt(x*x+y*y)*y);
}