11. 11章練習問題の答え・3

#include<stdio.h>
#include<math.h>
#define T 0.5     // calculation duration
#define max 10    // the maximum number of i
#define dt 0.01   // delta t
#define dx 0.1    // delta x
#define omega 1.2 // omega
#define a 1.0     // a
#define eps 0.00001
#define pi 3.141592
double  u[max+1][2]; //u_i^n is written as u[i][n], while n=0 or 1

int main(void){
  double uN,time=0;
  double err,r=a*dt/dx/dx; int i,j;
  for(i=0;i<=max;i++)
    u[i][0]=0.0;u[i][1]=0.0;
  for(i=0;i<=max;i++){     // define the initial condition
    u[i][0]=sin(pi*(double)i*dx);
  }
  u[0][1]=0; u[max][1]=0;  // put the boundary condition
  while (T>time){
   do {
     err=0.0;
     for(i=1;i<max;i++){
       uN=u[i][1];
       u[i][1]=(1-(double)omega)*uN
                + (double)omega*(1/(1+2*r)*u[i][0]
                + r/(1+2*r)*(u[i+1][1]+u[i-1][1]));
       if (err<fabs(u[i][1]-uN))
         err=fabs(u[i][1]-uN);
     }
   }while(err>eps);
   printf("%3.3f\t",time);
   time=time+dt;
   for(i=0;i<=max;i++)
     printf("%3.3f\t",u[i][0]);
   printf("\n");
   for(i=1;i < max;i++)
     u[i][0]=u[i][1];
  }
  return (0);
}