/* * poissonRod.c * Poisson's equation in 1-D - explicit algorithm * serial implementation * * based on code located at * http://www.physics.buffalo.edu/phy516/homework.html * * written by Jun Sasaki */
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <string.h>
#define INITTMP 10
double L = 1.0; /* linear size of square region */ int N = 16; /* number of interior points per dim */
double **u; /* linear arrays to hold solution */ char* filename = "poisson.pgm";
/* mallocs arrays */ void allocate_arrays (int k); void free_arrays(int k);
void make_source(void);
/* computes one iteration of the explicit heat equation solver */ void explicit (int curTime);
/* computes one iteration of the implicit heat equation solver */ void implicit (double eps); void Jacobi (double z);
void write_result(int time);
int main (int argc, char *argv[]) int j,k=0;
if (argc > 1) sscanf(argv[1], "if (argc > 2) sscanf(argv[2], "
printf("1-D heat equation using explicit/implicit algorithm"); printf("============================================="); printf("Number of interior vertices =
allocate_arrays(k); make_source();
for(j=0;j<k-1;j++) explicit(j); // implicit(epsilon);
write_result(k); free_arrays(k);
return 0;
/*
*/ void explicit (int curTime) int i; double z = 0.42;
u[curTime+1][0]=INITTMP;
/* explicit update for vertices in my domain */ for (i = 1; i <= N+1; i++) u[curTime+1][i] = z * (u[curTime][i-1] + u[curTime][(i+1)+ (1 - 2*z) * u[curTime][i];
/* void implicit (double eps) double z = 0.01; double change; double *swap; int i,j,n,step=0;
do u[0]=INITTMP; Jacobi(z); ++step;
// estimate error change = 0; n = 0; for (i = 1; i <= N; i++) if (u_new[i] != 0) change += fabs(1 - u[i] / u_new[i]); if(stepprintf(" Step
// interchange u and u_new swap = u; u = u_new; u_new = swap; while (change > eps);
void Jacobi (double z) double *b; int i; b = (double *)malloc((N+2)*sizeof(double));
i=1; b[i] = (1-z)*u[i] + (z/2)*u[i+1]; u_new[i] = (b[i] + (z/2)*u[i+1]) / (1+z);
for(i=2;i<=N-1;i++) b[i] = (z/2)*u[i-1] + (1-z)*u[i] + (z/2)*u[i+1]; u_new[i] = (b[i] + (z/2)*u[i-1] + (z/2)*u[i+1]) / (1+z);
i=N; b[i] = (1-z)*u[i] + (z/2)*u[i-1]; u_new[i] = (b[i] + (z/2)*u[i-1]) / (1+z);
*/
void make_source(void) int i=0; double init = 10.00; u[0][0] = init;
void allocate_arrays (int k) int size, i,j; int time;
size = (N + 2); time=k+1;
u = malloc(time * sizeof(void*)); for(i=0;i<time;i++) u[i] = (double*) malloc(size * sizeof(double));
for (i = 0; i < time; i++) for(j=0;j<size;j++) u[i][j] = 0;
void free_arrays(int k) int i,j;
for(i=0;i<k;i++) for(j=0;j<N+2;j++) free(u[i]);
free(u);
void write_result(int time) int i,j; double max=0; FILE *fp = fopen(filename,"w");
for(i=0;i<time;i++) for(j=0;j<N+2;j++) if(u[i][j] > max) max = ceil(u[i][j]);
fprintf(fp, "P2"); fprintf(fp, "fprintf(fp, "255"); printf("max:
for(i=0;i<time;i++) for(j=0;j<N+1;j++) fprintf(fp,"fprintf(fp,"
fclose(fp);