3.2.0.1 tanshin-euler.c

/* tanshin-euler.c -- 微分方程式の初期値問題を Euler 法で解く */

#include <stdio.h>
#include <math.h>

double omega2; // ω^2

int main(void)
{
  /* 開始時刻と終了時刻 */
  double a = 0.0, b = 1.0;
  /* 変数と関数の宣言 */
  int N, j;
  double t,x,y,h,x0,y0,fx(double, double, double), fy(double, double, double);
  /* 初期値 */
  x0 = 1.0; y0 = 0.0;
  /* 角振動数ω */
  double omega = 1.0;
  omega2 = omega * omega;
  /* 区間の分割数 N を入力してもらう */
  printf("N="); scanf("%d", &N);
  /* 小区間の幅 */
  h = (b-a) / N;
  /* 開始時刻と初期値のセット */
  t=a;
  x=x0; y=y0;
  printf("t=%f, (x,y)=(%f,%f)\n", t, x, y);
  /* Euler 法による計算 */
  for (j = 0; j < N; j++) {
    x += h * fx(t,x,y);
    y += h * fy(t,x,y);
    t += h;
    //printf("t=%f, (x,y)=(%f,%f)\n", t, x, y);
    printf("t=%f, (x,y)=(%f,%f) (%f,%f)\n", t, x, y,
           x0 * cos(omega *t) + y0 / omega * sin(omega * t),
           - x0 * omega * sin(omega * t) + y0 * cos(omega * t));
  }

  return 0;
}

/* 微分方程式 x'=fx(t,x,y) の右辺の関数 fx の定義 */
double fx(double t, double x, double y)
{
  return y;
}

/* 微分方程式 y'=fy(t,x,y) の右辺の関数 fy の定義 */
double fy(double t, double x, double y)
{
  return - omega2 * x;
}