3. 3章の答え・3

#include <stdio.h>
#include <math.h>
#define N 10
#define h 1.0
double f(double t,double x, double y);
double g(double t,double x, double y);
int main(void) {
  int i;
  double x=1000.0;
  double y=100.0;
  double t=0.0;
  double kx1,ky1,kx2,ky2,kx3,ky3,kx4,ky4;
  for (i = 0; i < N; i++) {
    kx1=h*f(t,x,y);ky1=h*g(t,x,y);
    kx2=h*f(t+0.5*h,x+kx1*0.5,y+ky1*0.5);ky2=h*g(t+0.5*h,x+kx1*0.5,y+ky1*0.5);
    kx3=h*f(t+0.5*h,x+kx2*0.5,y+ky2*0.5);ky3=h*g(t+0.5*h,x+kx2*0.5,y+ky2*0.5);
    kx4=h*f(t+h,x+kx3,y+ky3);ky4=h*g(t+h,x+kx3,y+ky3);
    x=x+(kx1+2.0*kx2+2.0*kx3+kx4)/6.0;
    y=y+(ky1+2.0*ky2+2.0*ky3+ky4)/6.0;
    t=t+h;
    printf("%4.4lf\t %4.4lf\t %e\n",x,y,t);
    }
}
double f(double t,double x, double y) {
  return (0.01*x-0.00001*x*y);
}
double g(double t,double x, double y) {
  return (-0.02*y+0.00001*x*y);
}