課題３−４のプログラム例

ホイン法

プログラム例

2変数の微分方程式にだけ適用できるシンプルなプログラム

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

/* 関数定義 */
double func_x(double x, double y);
double func_y(double x, double y, double m);
void heun(double x0, double y0, double m, double h, double t_max,
	  double (*fx)(double, double), double (*fy)(double, double, double));

/* メイン */
int main() {

  double x0=2.0, y0=0.5, m, t_max= 50.0, dt=0.001;

  printf("係数 mとdtを与えてください. --> ");
  scanf("%lf %lf", &m, &h);

  cpgopen("/xserv");
  cpgpap(8.0,0.3);
  cpgenv(0.,t_max,0.0,4.0,0,0);

  heun(x0, y0, m, dt, t_max, func_x, func_y);

}

/* ホイン法 */
void heun(double x0, double y0, double m, double h, double t_max,
	  double (*fx)(), double (*fy)()) {
  double x, y, t=0.0, k1[2], k2[2];

  x = x0; y = y0;

  for (t = 0.0 ; t < t_max ; t+=h) {

    x0 = x; y0 =y;  // 線を引くために前の値をx0,y0にとっておく

    /**********説明**************
     * (x,y)を使って、xとyの増分k1[0], k1[1]を計算し、
     * x+h*k1[0], y+h*k1[1]を使って、k2[0], k2[1]を計算
     * k1,k2をつかって、新しいx,yを求める
     * テキストでのxは、この問題での変数xではなく、時間tに相当することに注意
     ****************************/
    k1[0] = fx(x,y);
    k1[1] = fy(x,y,m);
    k2[0] = fx(x+h*k1[0], y+h*k1[1]);
    k2[1] = fy(x+h*k1[0], y+h*k1[1]);
    x = x + h/2.0 * ( k1[0] + k2[0] );
    y = y + h/2.0 * ( k1[1] + k2[1] );

    printf("t=%f \t x=%f \t y=%f \n", t, x, y);

    cpgsci(2);
    cpgmove(t, x0);
    cpgdraw(t,x);
    cpgsci(3);
    cpgmove(t, y0);
    cpgdraw(t,y);
  }
}

/* 関数の定義 */
double func_x(double x, double y){
  return x * (1.0 - y);
}
double func_y(double x, double y, double m) {
  return m * y * (x-1.0);
}

課題3-5 (発展課題)

微分方程式の数（未知数の数）や様々な方程式に対応できるようオイラー法、ルンゲクッタ法の計算関数を下のように作成し、それを利用してロトカ・ヴォルテラ方程式の計算をしてみる。

作成の方針は以下のとおりである。

作図等はメイン関数で行う．（もっと複雑な図になる場合は，作図用の関数を別途作成し，メイン関数から呼び出すようにする．）
euler(), rk4()を、多数の変数に対応できるようにする．
- 変数の数を関数euler()やrk4()に引数で与える (size_of_x)
- rk4()の内部の途中計算でつかうk1, k2, .. は配列にする．
- 方程式の定義 func() は変数を配列として受け取り，計算結果は引数を通じて返す．
係数mは、一括してparamという配列で、euler()やrk4()を経由して渡す。こうすることでeuler()やrk4()を修正しないで、新しい方程式に適用できる。

課題3-5のプログラム例

#include "cpgplot.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
/*
 * euler, rk4を変数の数に依存しない形に変えたプログラム
 * rk4を利用した例
 */
/* 関数定義 */
void euler(int n, double *xin, double *xout, double h, void (*f)(double *, double *, double *), double *params);
void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params);

/*
 * xin 入力、xout 出力、h 時間刻み、微分方程式の右辺
 * xin, xout は同じ大きさの配列とする
 */

void func(double *xin, double *xout, double *params){
  *xout = *xin * (1.0 - *(xin+1));
  *(xout+1) = *params * *(xin+1) * (*xin - 1.0);
}

/* メイン */
int main() {

  double x[2], x0[2] = {2.0,0.5}, t_max= 50.0, h=0.001;
  double t0=0., t=0.;
  double prev_x[2];  // 軌跡の最高地点調査用
  double params[1];
  printf("係数 mとdtを与えてください. --> ");

  scanf("%lf %lf", params, &h);

  cpgopen("/xserv");
  cpgpap(8.0,0.3);
  cpgenv(0.,t_max,0.0,5.0,0,0);
  cpgsci(2);
  cpgslw(3);

  x[0] = x0[0]; x[1] = x0[0]; t0=0.0;
  for (t = 0.0 ; t < t_max ; t+=h) {
    x0[0] = x[0]; x0[1] =x[1];  // 線を引くために前の値をx0,y0にとっておく

    // euler(2, x0, x, h, func, params);
    rk4(2, x0, x, h, func, params);

    // printf("t=%f \t x=%f \t y=%f \n", t, x[0], x[1]);

    cpgsci(2);       // 色の設定
    cpgmove(t0, x0[0]); // (t0,x0)にペンをアップで移動
    cpgdraw(t,x[0]);    // (t,x)までペンダウンで移動（線を書く）
    cpgsci(3);       // 色の設定
    cpgmove(t0, x0[1]); // (t0,y0)にペンアップで移動
    cpgdraw(t,x[1]);    // (t,y)までペンダウンで移動

    t0 = t;
  }
}

/************* ここ以下は固定 **************/
/* オイラー法 */
void euler(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal;

  funcVal = (double *) malloc(size_of_x * sizeof(double));

  f(x0, funcVal, params);  // 関数値の計算
  for(i=0; i<size_of_x; i++) {
    x[i] = x0[i] + h * funcVal[i]; //オイラー法
  }

}

/* ルンゲ・クッタ法 */
void rk4(int size_of_x, double *x0, double *x, double h,  void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal[4];

  for(i=0;i<4;i++) {
   funcVal[i] = (double *) malloc(size_of_x * sizeof(double));
  }
  f(x0, funcVal[0], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0*funcVal[0][i];
  f(x, funcVal[1], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0* funcVal[1][i];
  f(x, funcVal[2],params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h* funcVal[2][i];
  f(x, funcVal[3], params);

  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/6.0 * (funcVal[0][i]+2.0*funcVal[1][i]+2.0*funcVal[2][i]+funcVal[3][i]);

}

軌跡を描くプログラム例

ついでに、軌跡の最高点の位置を標準出力に書きだすことも行う。

#include "cpgplot.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
/*
 * euler, rk4を変数の数に依存しない形に変えたプログラム
 * rk4を利用した例
 */


/* 関数定義 */
void euler(int n, double *xin, double *xout, double h, void (*f)(double *, double *, double *), double *params);
void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params);

/*
 * xin 入力、xout 出力、h 時間刻み、微分方程式の右辺
 * xin, xout は同じ大きさの配列とする
 * 係数は配列paramで受け取る
 */

void func(double *xin, double *xout, double *params){
  *xout = *xin * (1.0 - *(xin+1));
  *(xout+1) = *params * *(xin+1) * (*xin - 1.0);
}

/* メイン */
int main() {

  double x[2], x0[2] = {2.0,0.5}, t_max= 50.0, h=0.001;
  double t0=0., t=0.;
  double prev_x[2];  // 軌跡の最高地点調査用
  double params[1];
  printf("係数 mとdtを与えてください. --> ");

  scanf("%lf %lf", params, &h);

  cpgopen("/xserv");
  cpgpap(5.0,1.0);
  cpgenv(0.,5.0,0.0,5.0,0,0);
  cpgsci(2);
  cpgslw(3);

  x[0] = x0[0]; x[1] = x0[0]; t0=0.0;

  cpgmove(x[0],x[1]);
  while (1) {
	  prev_x[0] = x0[0]; prev_x[1] = x0[1];          // ひとつ前の値を覚えておく
    x0[0] = x[0]; x0[1] =x[1];  // 線を引くために前の値をx0,y0にとっておく

    //euler(2, x0, x, h, func, params);
    rk4(2, x0, x, h, func, params);

    if(x0[1] > prev_x[1] && x[1] < x0[1])  // 最高地点の値出力
		printf("t=%f \t x=%f \t y=%f \n", t, x0[0], x0[1]);

    cpgdraw(x[0],x[1]);    // (x[0],x[1])までペンダウンで移動

    t += h;
  }

}

/************* ここ以下は固定 **************/
/* オイラー法 */
void euler(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal;

  funcVal = (double *) malloc(size_of_x * sizeof(double));

  f(x0, funcVal, params);  // 関数値の計算
  for(i=0; i<size_of_x; i++) {
    x[i] = x0[i] + h * funcVal[i]; //オイラー法
  }

}

/* ルンゲ・クッタ法 */
void rk4(int size_of_x, double *x0, double *x, double h,  void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal[4];

  for(i=0;i<4;i++) {
   funcVal[i] = (double *) malloc(size_of_x * sizeof(double));
  }
  f(x0, funcVal[0], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0*funcVal[0][i];
  f(x, funcVal[1], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0* funcVal[1][i];
  f(x, funcVal[2],params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h* funcVal[2][i];
  f(x, funcVal[3], params);

  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/6.0 * (funcVal[0][i]+2.0*funcVal[1][i]+2.0*funcVal[2][i]+funcVal[3][i]);

}

ソースの分割

大きなプログラムになると、すべてを１つのファイルに書き込むことは非効率的ですし、チームでの開発が困難です。

そこで、仕様を確定しておき、複数のファイルに分けてコードを書きます：

関数定義はヘッダーファイルに書く。
共通に利用するファイルを別に書く。

その例を、上記の軌跡描画のプログラムを例に示します。

ヘッダーファイル diff_methods.h

/* 微分方程式数値計算 関数定義 */
void euler(int n, double *xin, double *xout, double h, void (*f)(double *, double *, double *), double *params);
void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params);

微分方程式解法関数 diff_methods.c

/**** 微分方程式の数値計算関数 ****/
#include <stdlib.h>
/* オイラー法 */
void euler(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal;

  funcVal = (double *) malloc(size_of_x * sizeof(double));

  f(x0, funcVal, params);  // 関数値の計算
  for(i=0; i<size_of_x; i++) {
    x[i] = x0[i] + h * funcVal[i]; //オイラー法
  }

}

/* ルンゲ・クッタ法 */
void rk4(int size_of_x, double *x0, double *x, double h,  void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal[4];

  for(i=0;i<4;i++) {
   funcVal[i] = (double *) malloc(size_of_x * sizeof(double));
  }
  f(x0, funcVal[0], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0*funcVal[0][i];
  f(x, funcVal[1], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0* funcVal[1][i];
  f(x, funcVal[2],params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h* funcVal[2][i];
  f(x, funcVal[3], params);

  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/6.0 * (funcVal[0][i]+2.0*funcVal[1][i]+2.0*funcVal[2][i]+funcVal[3][i]);

}

ロトカ・ヴォルテラモデルの数値計算・描画 prac3_4adv.c

上に作成したヘッダーファイルをincludeします。

#include "cpgplot.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "diff_methods.h"
/*
 * euler, rk4を変数の数に依存しない形に変えたプログラム
 * rk4を利用した例
* /

double m; /* グローバル変数　*/


/*
 * xin 入力、xout 出力、h 時間刻み、微分方程式の右辺
 * xin, xout は同じ大きさの配列とする
 */

void func(double *xin, double *xout, double *params){
  *xout = *xin * (1.0 - *(xin+1));
  *(xout+1) = *params * *(xin+1) * (*xin - 1.0);
}

/* メイン */
int main() {

  double x[2], x0[2] = {2.0,0.5}, t_max= 50.0, h=0.001;
  double t0=0., t=0.;
  double prev_x[2];  // 軌跡の最高地点調査用
  double params[1];
  printf("係数 mとdtを与えてください. --> ");

  scanf("%lf %lf", params, &h);

  cpgopen("/xserv");
  cpgpap(5.0,1.0);
  cpgenv(0.,5.0,0.0,5.0,0,0);
  cpgsci(2);
  cpgslw(3);

  x[0] = x0[0]; x[1] = x0[0]; t0=0.0;

  cpgmove(x[0],x[1]);
  while (1) {
	  prev_x[0] = x0[0]; prev_x[1] = x0[1];          // ひとつ前の値を覚えておく
    x0[0] = x[0]; x0[1] =x[1];  // 線を引くために前の値をx0,y0にとっておく

    // euler(2, x0, x, h, func,params);
    rk4(2, x0, x, h, func, params);

    if(x0[1] > prev_x[1] && x[1] < x0[1])  // 最高地点の値出力
		printf("t=%f \t x=%f \t y=%f \n", t, x0[0], x0[1]);


    cpgdraw(x[0],x[1]);    // (x[0],x[1])までペンダウンで移動

    t += h;
  }

}

複数ファイルのコンパイル

コンパイルには、ソースファイルを並べます。

cc -o prac3-4adv prac3-4adv.c diff_methods.c -lm -lpgplot -lcpgplot -lX11

ヘッダーファイルが、別のディレクトリにあるときは"-I ディレクトリ名"のようなオプションをつけます。

また、diff_methods.cは再配置可能なオブジェクトファイルにコンパイルしておき、"diff_methods.o"を使うこともできます。たくさんのオプジェクトファイルをまとめてライブラリとしておけば、それをリンクするというかたちで使うことができます。

TOYOKI Hiroyasu
Updated: 2018-11-26

Home / Prog3Home