1 1階の常微分方程式

1.1 4次のルンゲ・クッタ法

前回の授業で学習した4次のルンゲ・クッタ法のプログラムの例を以下に示す。使用方法 は、以下の通りである。

• 初期条件$x_0$と$y_0$を聞いてくるので、それぞれをキーボードから入力する。
• 計算終了の$x$の値を聞いてくるので、それをキーボードから入力する。
• 計算ステップ数を聞いてくるので、それをキーボードから入力する。
• 他の微分方程式を計算したければ、関数funcの計算を変える。

=1

#include <stdio.h>
#include <math.h>
#define IMAX 100001
double func(double x, double y);
int mk_plot_data(char *a, double x[], double y[], int n);
int mk_graph(char *f);

/*================================================================*/
/*       main function                                            */
/*================================================================*/
int main(){
  double x[IMAX], y[IMAX];
  double final_x, h;
  double k1, k2, k3, k4;
  char temp;
  int ncal, i;

  /*-- set initial condition and cal range --*/

  printf("\ninitial value x0 = ");
  scanf("%lf%c", &x[0], &temp);

  printf("\ninitial value y0 = ");
  scanf("%lf%c", &y[0], &temp);

  printf("\nfinal x = ");
  scanf("%lf%c", &final_x, &temp);

  do{
    printf("\nNumber of calculation steps( <= %d ) = ",IMAX-1);
    scanf("%d%c", &ncal, &temp);
  }while(ncal > IMAX-1);        

  /* -- size of calculation step -- */

  h=(final_x-x[0])/ncal;

        
  /* -- 4th Runge Kutta Calculation -- */

  for(i=0; i < ncal; i++){
    k1=h*func(x[i],y[i]);
    k2=h*func(x[i]+h/2.0, y[i]+k1/2.0);
    k3=h*func(x[i]+h/2.0, y[i]+k2/2.0);
    k4=h*func(x[i]+h, y[i]+k3);

    x[i+1]=x[i]+h;
    y[i+1]=y[i]+1.0/6.0*(k1+2.0*k2+2.0*k3+k4);
  }

  /* -- make a graph -- */

  mk_plot_data("out.txt", x, y, ncal);
  mk_graph("out.txt");

  return 0;
}

/*================================================================*/
/*       define function                                          */
/*================================================================*/
double func(double x, double y){
  double dydx;

  dydx=sin(x)*cos(x)-y*cos(x);

  return(dydx);
}

/*==========================================================*/
/*   make a data file                                       */
/*==========================================================*/
int mk_plot_data(char *a, double x[], double y[], int n){
  int i;
  FILE *out;

  out = fopen(a, "w");

  for(i=0; i<=n; i++){
    fprintf(out, "%e\t%e\n", x[i], y[i]);
  }

  fclose(out);
  
  return 0;
}

/*==========================================================*/
/*   make a graph                                           */
/*==========================================================*/
int mk_graph(char *f){
  FILE *gp;
        
  gp = popen("gnuplot -persist","w");
        
  fprintf(gp, "reset\n");
  fprintf(gp, "set terminal postscript eps color\n");
  fprintf(gp, "set output \"graph.eps\"\n");
  fprintf(gp, "set grid\n");

  /* ----  plat graph -----*/

  fprintf(gp, "plot \"%s\" using 1:2 with line\n",f);
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "replot\n");
                                        
  pclose(gp);

  return 0;
}

1.2 プログラムの説明

このプログラムの大まかな構成は、

• main関数は、計算条件の設定と4次のルンゲ・クッタ法による計 算、グラフ作成のルーチンの呼び出しを行っている。
• 関数funcは、解くべき微分方程式 $dy/dx=f(x,y)$の$f(x,y)$の 計算を行っている。
• 関数mk_plot_dataは、計算結果のxとyの配 列の値を、ファイル(ハードディスク)に保管している。
• 関数mk_graphは、ファイルからグラフを作成している。

である。

各行の役割は、以下の通りである。

2行

数学関数を使っているので、math.hをインクルードしている。

3行

計算に使う配列の大きさIMAXを定義している。IMAXは最 大計算ステップ数に1以上を加えた値にしなくてはならない。

4-7行

プログラマーが作成した関数のプロトタイプ宣言。

11行

main関数の始まり。

12-16行

main関数で使う変数や配列の宣言。

21-28行

メッセージを書き出し、計算に必要な初期値と終値を取り込む。

30-33行

計算ステップ数の入力。ただし、配列の範囲を超えた場合、再度、入力を 促すようになっている。

38行

計算ステップのサイズ$h$の計算。

43-51行

4次のルンゲ・クッタ法の計算。

56行

計算結果をファイルに格納する関数(サブルーチン)の呼び出し。

57行

ファイルに書かれたデータをgnuplotでグラフにする関数(サブルーチン)の 呼び出し。

59行

main関数の戻り値。このプログラムでは、使われていない。

60行

main関数の終わり。

to 140mm

66行

関数funcの始まり。

• 関数名は、func
• 入力引数は、倍精度実数xと倍精度実数y
• 戻り値は、倍精度実数

67行

関数funcで使う変数宣言。

69行

解くべき微分方程式 $dy/dx=f(x,y)$の$f(x,y)$の計算。計算結果は、変数 dydxに格納。

71

関数funcの戻り値を、変数dydxの値として与えている。

72

関数funcの終わり。

to 140mm

78行

関数mk_plot_dataの始まり。

• 関数名は、mk_plot_data
• 入力引数は、文字列のポインターa、倍精度実数配列 xとy
• 戻り値は、整数

79-80行

関数mk_plot_dataで使う変数宣言。型名FILEは、 ファイルを使うときに必要である。

82行

ポインターaが示す文字列の名前で、書き込みモードでファイルを 開く。この後、ファイルは、outというポインターをつかってアク セスする。

84-86行

ファイルへデータの書き込み。

88行

ファイルを閉じる。

59行

mk_plot_data関数の戻り値。このプログラムでは、使われていない。

60行

mk_plot_data関数の終わり。

to 140mm

96行

関数mk_graphの始まり。

• 関数名は、mk_graph
• 入力引数は、文字列のポインターa。これは、データファ イル名を表す。
• 戻り値は、整数

97行

関数mk_graphで使う変数宣言。型名FILEは、 パイプを使うときに必要である。

99行

書き込みモードでファイルでパイプを開く。これは、ターミナルで、 「gnuplot -persist」と打ち込んだのと同じである。

101 -111行

gnuplotにダブルクォーテションで囲んだ命令を与えている。

113行

パイプを閉じている。

115行

mk_graph関数の戻り値。このプログラムでは、使われていない。

116行

mk_graph関数の終わり。

ホームページ: Yamamoto's laboratory
著者: 山本昌志