..................の間がプログラム本体です
program 15
..................
#include<stdio.h>
#define NMAX 4
/* void型関数サブプログラムの例 
  行列とベクトルを与えるとその積を返す
  (program 7 を関数を用いて書き換えたもの)*/

void product(double [][NMAX], double [], double []);

main()
{
  double a[NMAX][NMAX],x[NMAX],y[NMAX];
  int i,j;

  /*  行列に値をいれる(値自身は意味ないです) */
  a[0][0]=0.5;  a[0][1]=1.0;  a[1][0]=0.8;  a[1][1]=2.0;

/*  ベクトルに値をいれる */
  x[0]=1.0;  x[1]=2.0;

/* 積を求める関数を呼び出す */
  product(a,x,y);

/* 出力  */
  printf("matrix a:\n");
  for( i=0; i<2; ++i){
    for( j=0; j<2; ++j){
      printf("  %f",a[i][j]);
    }
    printf("\n");
  }
  printf("vector x:\n");
  for( i=0; i<2; ++i){ 
    printf("  %f\n",x[i]);
  }
  printf("vector y=a*x:\n");
  for( i=0; i<2; ++i){ 
    printf("  %f\n",y[i]);
  }
}

void product(double m[][NMAX], double v[], double v2[])
{
 int i,j;

/*  ベクトルv2を初期化(値を0に) */
  for( i=0; i<2; ++i){
     v2[i]=0;
  }    

/*  m*vによりv2を求める */
  for( i=0; i<2; ++i){
    for( j=0; j<2; ++j){
       v2[i] += m[i][j]*v[j];
    }
  }
}
..................
program 16a
..................
#include<stdio.h>
#include<math.h>
/* 
void型関数の応用例
最低次のsimplectic integrator (対称分解)によって
 一次元調和振動子の運動の初期値問題を解く。
*/

void simplectic(double *, double *, double);

main()
{
 double delta_t,q,p;
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&q,&p);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  simplectic(&q,&p,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,q,p,(p*p+q*q)/2);
 }
}

/* 1ステップ分の計算をする関数 */
void simplectic(double *q, double *p, double delta_t)
{  
  *q += *p * delta_t/2;
  *p -= *q * delta_t;
  *q += *p * delta_t/2;
  }
..................
program 16b
..................
#include<stdio.h>
#include<math.h>
/* 
16aと同じ問題を、配列を関数に渡す方式で書き直すとこうなる
*/

void simplectic(double *, double);

main()
{
 double delta_t,qp[2];
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&qp[0],&qp[1]);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  simplectic(qp,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,qp[0],qp[1],
 (qp[0]*qp[0]+qp[1]*qp[1])/2);
 }
}

/* 1ステップ分の計算をする関数 */
void simplectic(double *qp, double delta_t)
{  
  qp[0] += qp[1] * delta_t/2;
  qp[1] -= qp[0] * delta_t;
  qp[0] += qp[1] * delta_t/2;
  }
..................
program 16c
..................
#include<stdio.h>
#include<math.h>
/* 
一見、こちらのプログラムを書きたくなるが、
Cの場合、このやりかたではうまくいかない(16aとの違いに注目)
*/

void simplectic(double , double , double);

main()
{
 double delta_t,q,p;
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&q,&p);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  simplectic(q,p,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,q,p,(p*p+q*q)/2);
 }
}

/* 1ステップ分の計算をする関数 */
void simplectic(double q, double p, double delta_t)
{  
  q += p * delta_t/2;
  p -= q * delta_t;
  q += p * delta_t/2;
  }
..................
program 17
..................
#include<stdio.h>
#include<math.h>
/* 
16より少しだけ汎用のプログラム
*/

void simplectic(double *, double *, double);
double hq(double, double);
double hp(double, double);

main()
{
 double delta_t,q,p;
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&q,&p);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  simplectic(&q,&p,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,q,p,(p*p+q*q)/2);
 }
}

/* 1ステップ分の計算をする関数 */
void simplectic(double *q, double *p, double delta_t)
{  
  *q += hq(*q,*p) * delta_t/2;
  *p += hp(*q,*p) * delta_t;
  *q += hq(*q,*p) * delta_t/2;
  }

/* simplectic integratorの座標部分
 (調和振動子の場合はほとんど自明 */
double hq(double q, double p)
{
 return p;
}
/* simplectic integratorの運動量部分
 (調和振動子の場合はほとんど自明 */
double hp(double q, double p)
{
 return -q;
}
..................
program 18
..................
#include<stdio.h>
#include<math.h>
/* 
16と同じ問題を四次のRunge-Kuttaで解くプログラム
*/

void rungekutta(double *, double *, double);
void onestep(double, double, double *, double *, double);

main()
{
 double delta_t,q,p;
 int i,step;

/* 時間刻み、ステップ数、初期の座標と運動量を入力 */
 scanf("%lf%d%lf%lf",&delta_t,&step,&q,&p);

/* 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ */
 for(i=1; i<step; ++i){
  rungekutta(&q,&p,delta_t);

/* 各ステップごとに、時刻、座標、運動量、エネルギーを出力 */
  printf("%f  %f   %f   %f\n",i*delta_t,q,p,(p*p+q*q)/2);
 }
}

/* 1ステップ分の計算をする関数 */
void rungekutta(double *x, double *v, double delta_t)
{ 
  double k1x,k1v,k2x,k2v,k3x,k3v,k4x,k4v;
  onestep(*x, *v, &k1x, &k1v,delta_t);
  onestep(*x+k1x/2, *v+k1v/2, &k2x, &k2v, delta_t);
  onestep(*x+k2x/2, *v+k2v/2, &k3x, &k3v, delta_t);
  onestep(*x+k3x, *v+k3v, &k4x, &k4v, delta_t);
  *x += k1x/6 + k2x/3 + k3x/3 + k4x/6;
  *v += k1v/6 + k2v/3 + k3v/3 + k4v/6;
  }

/* rungekuttaの一段分
 (調和振動子の場合はほとんど自明 */
void onestep(double x, double v, double *xx, double *vv, double delta_t)
{
 *xx = v*delta_t;
 *vv = -x*delta_t;
}