気になる 番組やタレントやCMやモデルをピックアップ

/**********************************************************************
y' = x+ y, 0<=x<=1,y(0)=1の常微分方程式の解をもとめる
 **********************************************************************/
#include <stdio.h>
#define N 20
void runge_kutta (double x,double y,double a,double b,int n,double (*f)(double,double));
double f1(double x,double y);

int main(void)
{
  double y0=1.0,a=0.0,b=1.0;
  runge_kutta(0.0, y0, a, b, N, f1);
  return 0;
}

void runge_kutta (double x,double y,double a,double b,int n,double (*f)(double,double))
{
  int i;
  double h,k1,k2,k3,k4;
  h = (b-a)/n;
  x = a;
  printf("x=%lf y=%lf\n",x,y);/* 初期値を出力 */
  for(i=0;i<=n-1;i++){
    k1 = f(x,y);
    k2 = f(x + h/2, y + h/2 * k1);
    k3 = f(x + h/2, y + h/2 * k2);
    k4 = f(x + h, y + h* k3);
    y += h/6 * (k1 + 2*k2 + 2*k3 + k4);
    x += h;
    printf("x=%lf y=%lf\n",x,y);
  }
}

double f1(double x,double y)
{
  return x+y;
}
/**********************************************************************
$ ./a.out
x=0.000000 y=1.000000
x=0.050000 y=1.052542
x=0.100000 y=1.110342
x=0.150000 y=1.173668
x=0.200000 y=1.242805
x=0.250000 y=1.318051
x=0.300000 y=1.399718
x=0.350000 y=1.488135
x=0.400000 y=1.583649
x=0.450000 y=1.686624
x=0.500000 y=1.797442
x=0.550000 y=1.916506
x=0.600000 y=2.044237
x=0.650000 y=2.181082
x=0.700000 y=2.327505
x=0.750000 y=2.484000
x=0.800000 y=2.651082
x=0.850000 y=2.829294
x=0.900000 y=3.019206
x=0.950000 y=3.221419
x=1.000000 y=3.436563
**********************************************************************/

/**********************************************************************
y' = x+ y, 0<=x<=1,y(0)=1の常微分方程式の解をもとめる
 **********************************************************************/
#include <stdio.h>
#define N 20

double f1(double x,double y);
void heun(double x,double y,double a,double b, int n, double(*f)(double,double));

int main(void)
{
  heun(0.0, 1.0, 0.0, 1.0, N, f1 );
  return 0;
}
void heun(double x,double y,double a,double b,int n,double (*f)(double,double))
{
  int i;
  double h,k1,k2;
  h = (b-a)/n;
  x = a;
  for(i =0; i<=n-1; i++){
    k1 = f(x,y);
    k2 = f(x + h, y + h*k1 );
    y += h / 2 * ( k1 + k2);
    x += h;
    printf("x=%lf y=%lf\n",x,y);
  }
}

double f1(double x,double y)
{
  return x+y;
}

/**********************************************************************
$ ./a.out
x=0.050000 y=1.052500
x=0.100000 y=1.110253
x=0.150000 y=1.173529
x=0.200000 y=1.242609
x=0.250000 y=1.317793
x=0.300000 y=1.399393
x=0.350000 y=1.487736
x=0.400000 y=1.583170
x=0.450000 y=1.686058
x=0.500000 y=1.796781
x=0.550000 y=1.915741
x=0.600000 y=2.043360
x=0.650000 y=2.180082
x=0.700000 y=2.326374
x=0.750000 y=2.482726
x=0.800000 y=2.649653
x=0.850000 y=2.827698
x=0.900000 y=3.017430
x=0.950000 y=3.219448
x=1.000000 y=3.434382
 **********************************************************************/