Want more features on Pastebin? Sign Up, it's FREE!
Guest

jetxee

By: a guest on Jan 20th, 2009  |  syntax: C  |  size: 2.14 KB  |  views: 135  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
This paste has a previous version, view the difference. Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5.  
  6. // overwrite u on success
  7. double* stepEuler(double* u, int const n, double const mu) {
  8.         int i;
  9.         double *un=malloc(sizeof(double)*n);
  10.         if (!un) {
  11.                 return 0;
  12.         }
  13.         if (!u || n < 3) {
  14.                 free(un);
  15.                 return 0;
  16.         }
  17.         for (i=1; i<(n-1); ++i) { // inner points
  18.                 un[i]=u[i]+mu*(u[i-1]+u[i+1]-2*u[i]);
  19.         }
  20.         un[0]=u[0]+mu*2*(u[1]-u[0]); // zero-flux on boundaries
  21.         un[n-1]=u[n-1]+mu*2*(u[n-2]-u[n-1]);
  22.         memcpy(u,un,sizeof(double)*n);
  23.         free(un);
  24.         return u;
  25. }
  26.  
  27. double* makeu0(int const n) {
  28.         double *u0=0;
  29.         u0=malloc(sizeof(double) * n);
  30.         if (u0) {
  31.                 int i;
  32.                 for (i=0; i<n; ++i) {
  33.                         double xi=M_PI*i/n;
  34.                         u0[i]=sin(xi)*sin(xi);
  35.                 }
  36.         }
  37.         return u0;
  38. }
  39.  
  40. double sum_u(double const *u, int const n) {
  41.         if ((u) && (n > 1)) {
  42.                 double sum=0.0;
  43.                 int i;
  44.                 for (i=0; i<n; ++i) {
  45.                         sum += u[i];
  46.                 }
  47.                 return sum;
  48.         } else {
  49.                 return 0.0;
  50.         }
  51. }
  52.  
  53. void print_u(double const *u, int const n) {
  54.         if ((u) && (n > 1)) {
  55.                 int i;
  56.                 for (i=0; i<n; ++i) {
  57.                         printf("%g,",u[i]);
  58.                 }
  59.                 printf("\n");
  60.         }
  61. }
  62.  
  63. int main(int argc, char *argv[]) {
  64.         int n;
  65.         double *u0, *ret;
  66.         if (argc != 2) {
  67.                 printf("usage: ceuler n>=3\n");
  68.                 return 1;
  69.         } else {
  70.                 n=atoi(argv[1]);
  71.                 if (n < 3) {
  72.                         return 0;
  73.                 }
  74.         }
  75.         u0=makeu0(n);
  76.         ret=stepEuler(u0,n,0.5);
  77.         if (ret) {
  78.                 printf("%g\n",sum_u(u0,n));
  79.                 free(u0);
  80.         } else {
  81.                 printf("failure\n");
  82.                 free(u0);
  83.                 return 2;
  84.         }
  85. }
clone this paste RAW Paste Data