/*  prog2b.c -- Solve viscous Burgers Equation
    Written Apr 6, 2003 by Eric Olson for Mathematics 484
 
    For information about the central finite difference scheme used
    here see chapter 7 in Numerical Analysis for Applied Mathematics
    Science and Engineering by Greenspan and Casulli.  Information 
    about the heat and Burgers equations may be found in 13 and 14 in 
    A First Course in the Numerical Analysis of Differential Equations 
    by Iserles. */

#define N 16
#include <math.h>
#include <stdio.h>

double u[N],f[N],dx;

init(double u[N]){
    int i;
    for(i=0;i<N;i++){
        double x=dx*i;
        double xx=M_PI-x; xx*=xx;
        u[i]=8*exp(-xx);
    }
}

force(double f[N],double u[N]){
    int i;
    for(i=0;i<N;i++) {
        int ip1=(i+1)%N; 
        int im1=(i+N-1)%N;
        f[i]=(-u[i]*(u[ip1]-u[im1])+(u[im1]-2*u[i]+u[ip1])/dx)/dx;
    }
    return;
}

euler(double u[N],double tn,int n){
    double dt=tn/n;
    int i,j;
    for(j=0;j<n;j++) {
        force(f,u);
        for(i=0;i<N;i++) {
            u[i]=u[i]+dt*f[i];
        }
    }
    return;
}

vectwr(double u[N],double t){
    int i;
    printf("\n\n#At time %g...\n",t);
    for(i=0;i<N;i++){
        printf("%g %g\n",dx*i,u[i]);
    }
}

main(){
    double t;
    int i;
    dx=2*M_PI/N;
    printf("#prog2b.c -- Compute viscous Burger's Equation\n");
    
    init(u);
    vectwr(u,0.0);

    t=0.0;
    for(i=1;i<=8;i++){
        euler(u,0.125,4096);
        vectwr(u,i*0.125);
    }
    return 0;
}