Sunday, February 17, 2013

C code - Poisson Equation by finite difference method



//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");
            }
}

5 comments:

  1. Someone asked explanation for this code. I am putting them with his question here for your benefit.
    some 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.

    ReplyDelete
    Replies
    1. how to do this if right hand is dirac delta function

      Delete
    2. sir,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.

      Delete
  2. how g can be changed to a dirac delta function

    ReplyDelete
  3. Hello Sir,
    can 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

    ReplyDelete