//program to
solve Poisson's equation uxx+uyy=g(x,y)

// by Mahesha M G, MIT

#include
<stdio.h>

#define n
5 //dimension rows

#define m
4 //dimension columns

#define
gx(x,y) -4

main()

{

int i,j,k;

float
a[10][10],g[10][10],h,xl,xh,yl,x;

FILE *fp;

clrscr();

fp=fopen("c:\\poisson.dat","w");

printf("enter
xl: ");

scanf("%f",&xl);

printf("enter
xh: ");

scanf("%f",&xh);

printf("enter
yl: ");

scanf("%f",&yl);

h=(xh-xl)/(n-1);

for(i=1;i<=n;i++)

{

x=xl;

for(j=1;j<=m;j++)

{

g[i][j]=gx(x,yl);

x+=h;

}

yl+=h;

}

printf("enter
boundary conditions\n");

for(i=1;i<=n;i++)

{

for(j=1;j<=m;j++)

{

if(i==1||i==n||j==1||j==m)

{

printf("a[%d,%d] = ",i,j);

scanf("%f",&a[i][j]);

}

}

}

for(i=2;i<n;i++)

for(j=2;j<m;j++)

a[i][j]=a[n][j];
//initialization ends

for(k=0;k<100;k++)

{

for(i=2;i<n;i++)

{

for(j=2;j<m;j++)

{

a[i][j]=(a[i-1][j]+a[i+1][j]+a[i][j-1]+a[i][j+1]-(h*h*g[i][j]))/4;

}

}

}
//calculation by
Gauss-Seidel Method

for(i=1;i<=n;i++)

{

for(j=1;j<=m;j++)

fprintf(fp,"%f\t",a[i][j]);

fprintf(fp,"\n");

}

}

Someone asked explanation for this code. I am putting them with his question here for your benefit.

ReplyDeletesome doubts regarding this code..

1. What is the purpose of the variables xl,xh,yl?? What should we give input to them?? Earlier i thought that we are supposed to give the lower and higher limits to "x" and then we are calculating step size "h". but then i noted we are specifying the dimensions of rows and columns in "m" and "n". Then what are we specifying in xl,xh,yl??

2.In poisson's equation, we have a charge distribution rho which is given and by solving poisson we can tell the potential. In this code i think we are specifying the rho with gx(x,y) = - 4. So even if we change the value of "-4" to "0" there is no change in output. So i am having trouble to understand this.

ANSWER:

I agree that the explanation part is missing in the code.

Q1: Xl is lower value of x and xh is higher value of x. (xh-xl)/(n-1) gives step size. here n is number of grid points along the row. Similarly yl is the lower value of y. Upper value will be decided by code. I made step size along y same as that along x so that the finite difference expression becomes simple.

Q2: Let us take the following example:

Torsion on a rectangular bar subject to twisting is governed by del 2T = - 4. Given the condition T = 0 on boundary, find T over a cross section of a bar of size 9 cm X 12 cm. Use a grid size of 3 cm X 3 cm.

Inputs:

xl=0

xh=9

yl=0

all boundary points a(i,j) = 0

Following is the output.

0.000000 0.000000 0.000000 0.000000

0.000000 11.571428 11.571428 0.000000

0.000000 14.464285 14.464285 0.000000

0.000000 11.571428 11.571428 0.000000

0.000000 0.000000 0.000000 0.000000

If you make g(x,y)=0, then it reduces to Laplace equation and for same input conditions, output is following.

0.000000 0.000000 0.000000 0.000000

0.000000 0.000000 0.000000 0.000000

0.000000 0.000000 0.000000 0.000000

0.000000 0.000000 0.000000 0.000000

0.000000 0.000000 0.000000 0.000000

So there is difference. I wrote the code to solve the above problem which is taken from Numerical methods by Balagurusamy. We can change g(x,y). Also you give step size instead of defining the dimension of matrix as mXn.

Mail me if you have any further query.

how to do this if right hand is dirac delta function

Deletesir,i need the c source code for blasius solution for flat plate using finite difference method would u plz give me because i m from mechanical m.tech final year doind research work in iit. so kindly send it to my email address ranjan333999@gmail.com sir i request you plz kindly do it as soon as possible.

Deletehow g can be changed to a dirac delta function

ReplyDeleteHello Sir,

ReplyDeletecan you please give me any idea about how to solve elelctrostatic poisson equation having periodic boundary condition in the x direction

and fixed(Dirichlet) boundary condition in the y direction..?

Thanks