jetxee

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// overwrite u on success
double* stepEuler(double* u, int const n, double const mu) {
        int i;
        double *un=malloc(sizeof(double)*n);
        if (!un) {
                return 0;
        }
        if (!u || n < 3) {
                free(un);
                return 0;
        }
        for (i=1; i<(n-1); ++i) { // inner points
                un[i]=u[i]+mu*(u[i-1]+u[i+1]-2*u[i]);
        }
        un[0]=u[0]+mu*2*(u[1]-u[0]); // zero-flux on boundaries
        un[n-1]=u[n-1]+mu*2*(u[n-2]-u[n-1]);
        memcpy(u,un,sizeof(double)*n);
        free(un);
        return u;
}

double* makeu0(int const n) {
        double *u0=0;
        u0=malloc(sizeof(double) * n);
        if (u0) {
                int i;
                for (i=0; i<n; ++i) {
                        double xi=M_PI*i/n;
                        u0[i]=sin(xi)*sin(xi);
                }
        }
        return u0;
}

double sum_u(double const *u, int const n) {
        if ((u) && (n > 1)) {
                double sum=0.0;
                int i;
                for (i=0; i<n; ++i) {
                        sum += u[i];
                }
                return sum;
        } else {
                return 0.0;
        }
}

void print_u(double const *u, int const n) {
        if ((u) && (n > 1)) {
                int i;
                for (i=0; i<n; ++i) {
                        printf("%g,",u[i]);
                }
                printf("\n");
        }
}

int main(int argc, char *argv[]) {
        int n;
        double *u0, *ret;
        if (argc != 2) {
                printf("usage: ceuler n>=3\n");
                return 1;
        } else {
                n=atoi(argv[1]);
                if (n < 3) {
                        return 0;
                }
        }
        u0=makeu0(n);
        ret=stepEuler(u0,n,0.5);
        if (ret) {
                printf("%g\n",sum_u(u0,n));
                free(u0);
        } else {
                printf("failure\n");
                free(u0);
                return 2;
        }
}