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

ホイン法

方程式の右辺をそれぞれ別の関数として定義した標準例

#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("/CGW");
  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);
}

発展課題として説明したプログラム例

次の課題を行う上で，このようなプログラムを理解しましょう．

作図等はメイン関数で行う．（もっと複雑な図になる場合は，作図用の関数を別途作成し，メイン関数から呼び出すようにする．）
euler(), rk4()を、多数の変数に対応できるようにする．
- 変数の数を関数euler()やrk4()に引数で与える (size_of_x)
- rk4()の内部の途中計算でつかうk1, k2, .. は配列にする．
- 方程式の定義 func() は変数を配列として受け取り，計算結果は引数を通じて返す．
係数mはeuler()やrk4()を経由しないように関数funcに与える．そうすることにより，別のモデル方程式を扱うとき，rk4()などを変更しないですむ．
- 下記では，mをグローバル変数とした．

横軸をtとした描画（3-2の拡張)

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

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

/* 関数定義 */
void euler(int n, double *xin, double *xout, double h, void (*f)(double *, double *));
void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *));
/*
 * xin 入力、xout 出力、h 時間刻み、微分方程式の右辺
 * xin, xout は同じ大きさの配列とする
 */

/**
  ロトカ・ヴォルテラ微分方程式
*/
void func(double *xin, double *xout){
  *xout = *xin * (1.0 - *(xin+1));
  *(xout+1) = m * *(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.;

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

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

  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);
    rk4(2, x0, x, h, func);

    // 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;
  }

}

/************* ここ以下は触らない **************/
/**
  【変数の数によらずに使える微分方程式の数値解法関数の定義】
  ・微分方程式の変数を配列として、配列のポインタx0を与え、新しい値をx1に代入してもらう
  ・変数の数は size_of_x に与える
 **/

/* オイラー法 */
void euler(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *)) {

  int i;
  double *funcVal;

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

  f(x0, funcVal);  // 関数値の計算
  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 *)) {

  int i;
  /* double funcVal[4][size_of_x]; */ /* VC++ではエラーになる->mallocで確保 */

  double *funcVal[4]; /* 配列確保 */
  for(i=0;i<4;i++) {
   funcVal[i] = (double *) malloc(size_of_x * sizeof(double)); //変数の数だけのメモリを確保
  }

  f(x0, funcVal[0]); // 現在の値の配列x0を渡して，微分方程式の右辺を計算しfuncValに代入してもらう
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0*funcVal[0][i];
  f(x, funcVal[1]);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0*funcVal[1][i];
  f(x, funcVal[2]);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h*funcVal[2][i];
  f(x, funcVal[3]);

  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を利用した例
*; /

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

/* 関数定義 */
void euler(int n, double *xin, double *xout, double h, void (*f)(double *, double *));
void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *));
/*
 * xin 入力、xout 出力、h 時間刻み、微分方程式の右辺
 * xin, xout は同じ大きさの配列とする
 */

void func(double *xin, double *xout){
  *xout = *xin * (1.0 - *(xin+1));
  *(xout+1) = m * *(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];  // 軌跡の最高地点調査用

  printf("係数 mとdtを与えてください. --> ");

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

  cpgopen("/CGW");
  cpgpap(5.0,1.0);
  cpgenv(0.,5.0,0.0,5.0,0,0);

  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);
    rk4(2, x0, x, h, func);

    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 *)) {

  int i;
  double *funcVal;

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

  f(x0, funcVal);  // 関数値の計算
  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 *)) {

  int i;
  double *funcVal[4];

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

  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]);

}

TOYOKI Hiroyasu
Updated: 2016-12-07

Home / Prog3Home