2階常微分方程式

$z=y'=\frac{{\rm d}}{{\rm d}x}y$と置くのでこの部分はすべて共通。

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

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

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;                  /* 刻み幅 */

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

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

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