文書の過去の版を表示しています。

2階常微分方程式

ダウンロード版

オイラー法-2階版

#include <stdio.h>
#include <math.h>
double dydx(double x, double y, double z);
double dzdx(double x, double y, double z);
int main(void) {
  int i, n = 100;                    /* 繰り返し回数 */
  double x = 0.0, y = 0.0, z = 0.0;  /* 初期値 */
  double dx = 0.01;                  /* 刻み幅 */
  double dy, dz;
  printf("%f %f %f\n", x, y, z);
  for (i=0; i<n; i++) {
    dy = dydx(x, y, z)*dx;
    dz = dzdx(x, y, z)*dx;
    x = x + dx;
    y = y + dy;
    z = z + dz;
    printf("%f %f %f\n", x, y, z);
  }
  return 0;
}

修正オイラー法-2階版

#include <stdio.h>
#include <math.h>
double dydx(double x, double y, double z);
double dzdx(double x, double y, double z);
int main(void) {
  int i, n = 100;                   /* 繰り返し回数 */
  double x = 0.0, y = 0.0, z = 0.0; /* 初期値 */
  double dx = 0.01;                 /* 刻み幅 */
  double dy1, dy2, dz1, dz2;
  printf("%f %f %f\n", x, y, z);
  for (i=0; i<n; i++) {
    dy1 = dydx(x, y, z)*dx;
    dz1 = dzdx(x, y, z)*dx;
    dy2 = dydx(x + dx, y + dy1, z +dz1)*dx;
    dz2 = dzdx(x + dx, y + dy1, z +dz1)*dx;
    x = x + dx;
    y = y + (dy1 + dy2)/2.0;
    z = z + (dz1 + dz2)/2.0;
    printf("%f %f %f\n", x, y, z);
  }
  return 0;
}

ルンゲ・クッタ法-2階版

#include <stdio.h>
#include <math.h>
double dydx(double x, double y, double z);
double dzdx(double x, double y, double z);
int main(void) {
  int i, n = 100;                    /* 繰り返し回数 */
  double x = 0.0, y = 0.0, z = 0.0;  /* 初期値 */
  double dx = 0.01;                  /* 刻み幅 */
  double dy1, dy2, dy3, dy4, dz1, dz2, dz3, dz4; 
  printf("%f %f %f\n", x, y, z);
  for (i=0; i<n; i++) {
    dy1 = dydx(x, y, z)*dx;
    dz1 = dzdx(x, y, z)*dx;
    dy2 = dydx(x + dx/2., y + dy1/2., z +dz1/2.)*dx;
    dz2 = dzdx(x + dx/2., y + dy1/2., z +dz1/2.)*dx;
    dy3 = dydx(x + dx/2., y + dy2/2., z + dz2/2.)*dx;
    dz3 = dzdx(x + dx/2., y + dy2/2., z + dz2/2.)*dx;
    dy4 = dydx(x + dx, y + dy3, z + dz3)*dx;
    dz4 = dzdx(x + dx, y + dy3, z + dz3)*dx;
    x = x + dx;
    y = y + (dy1 + 2.0*(dy2 + dy3) + dy4)/6.0;
    z = z + (dz1 + 2.0*(dz2 + dz3) + dz4)/6.0;
    printf("%f %f %f\n", x, y, z);
  }
  return 0;
}

共通部分

<latex>z=y'=\fracrm_d{{\rm d}x}y</latex>と置くのでこの部分はすべて共通。

double dydx(double x, double y, double z) {
  return z;
}

問題によって変更する部分

<latex>x, y, y'(=z), y</latex>で書かれた式を<latex>y=\cdots</latex>の形に直して、その右辺の式を計算する関数をここに定義する。例は、<latex>y+2y'+y=1</latex>すなわち<latex>y=1-2y'-y=1-2z-y</latex>の場合の例になっている。

double dzdx(double x, double y, double z) {
  double ret;
  ret = 1.0 - 2.0*z - y;
  return ret;
}

また、問題や用途によって「繰り返し回数」、「初期値」、「刻み幅」は適当なものに適宜変更せよ。

  int i, n = 100;                    /* 繰り返し回数 */
  double x = 0.0, y = 0.0, z = 0.0;  /* 初期値 */
  double dx = 0.01;                  /* 刻み幅 */

@surface

目次