Solving Poisson's Equation Using 5-Point Stencil In A Rectangular Domain

Solving Poisson's Equation Using 5-Point Stencil:

Example:Consider the Poisson's equation on a unit square domain with the source term given by $q(x)=100*\sin(\pi x)*\sin(\pi y)$. The boundaries of the domain are maintained at $u=0$. Reactangular region: $0 \leq x,y \leq 1$. Determine the solution using finite difference method.

Solution:5-point stencil and Gauss-Seidel methods are used to solve Poisso's equations. Here a Mathematica code is give to sove this following program. dx = 0.1;
 x = Range[0, 1, dx];
 y = x;
 n = Length[x];m = Length[y];
 xy = Transpose[Join[{x}, {y}]];
(*Print["xy: ",xy//MatrixForm];*)
 xMeshgrid = ConstantArray[10, {n, n}];
 yMeshgrid = ConstantArray[10, {m, m}];
 o[Do[xMeshgrid[[i, j]] = x[[i]];
 yMeshgrid[[j, i]] = y[[i]], {j, 1, m}], {i, 1, n}];
 (*Print["xMeshgrid: ",xMeshgrid//MatrixForm];
 Print["yMeshgrid: ",yMeshgrid//MatrixForm];*)
 X = Transpose[xMeshgrid];
 Y = Transpose[yMeshgrid];
 (*Print["X: ",X//MatrixForm];
Print["Y: ",Y//MatrixForm];*)
 XY = Transpose[Join[{Flatten[xMeshgrid]}, {Flatten[yMeshgrid]}]];
 (*Print["meshgrid XY: ",XY//MatrixForm];*)
 h = 1/(n - 1); k = 1/(m - 1);
 q = 100*Sin[Pi X]*Sin[Pi Y];
 q[[-1, ;;]] = 0;
 q[[;; , -1]] = 0;
 Print["q : ", q // MatrixForm];
 u = ConstantArray[0, {n, n}];
u0 = u;
res = 1;
tol = 10^-6;
iter = 0;

 While[res > tol, iter = iter + 1;(*Print["Iteration: ",iter];*)
  Do[Do[u[[i, j]] = ((u[[i + 1, j]] + u[[i - 1, j]])/h^2 + (u[[i, j + 1]] + u[[i, j - 1]])/k^2 + q[[i, j]])/(2/h^2 + 2/k^2),
    {j, 2, n - 1}],
   {i, 2, n - 1}]; res = Max[Abs[u0 - u]];
 (*Print["res: ",res];*)
 u0 = u];
 Print["Total iteration requird for [", n, "x", n, "] grid is: ", iter];
 Print["solution u: ", u0 // MatrixForm]
 data = MapThread[Append, {XY, Flatten[u0]}]
 p1 = ListContourPlot[data, ContourLabels -> All, FrameLabel -> {"x", "y"}, PlotLegends -> Automatic]

Here contour plot of this solution is attached below:

  • By changing the value of dx, you can change the grid size. Warning! If you increase grid size from 21by21 to 81by81 or more it will take very large number of iterations to converge the solution. Here Gauss-Seidel method is applied.
  • Mathematica Official website: Visit Mathematica
  • Mathematica Documentation Center: Visit Mathematica Documentation Center
Location: India

Comments

Post a Comment

If you have any doubts, please let me know!

Popular posts from this blog

Solve Poisson's Equation in Matlab

Image
Solution Poisson's Equation In a Rectangular Domain Using 5-Point Stencil Poisson's equation in a square domain $0 \leq x,y \leq 1$ with source term given by $q(x)=100*\sin(\pi x)*\sin(\pi y)$. The boundaries of the domain are maintained at $u=0$. Determine the solution using finite difference method. Solution: Discretizing domain uniformly with step size h and k along x and y direction. Using five point stencil, equation at each node is given by $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}+\frac{u_{i,j+1}-2u{i,j}+u{i,j-1}}{k^{2}}=-q(i,j)$$ $$u_{i,j}=\frac{1}{\frac{2}{h^{2}}+\frac{2}{k^{2}}}(\frac{u_{i+1,j}+u_{i-1,j}}{h^{2}}+\frac{u_{i,j+1}+u{i,j-1}}{k^{2}}+q(i,j))$$ where $q(i,j)=100*\sin(\pi x_{i})*\sin(\pi y_{j})$ Here a MATLAB code is given below. Choosen $m=n=81$, these are the grid spacing. cme_assignment8_b clear; close; clc; n = 81; m = 81; x = 0:1/(n-1):1; y = 0:1/(m-1):1; h = 1/(n-1); k = 1/(m-1); [Y,X] = meshgrid(y,x); q = 100*...
Read more

Backward Substitution in MATLAB

Backward Substitution Backward Substitution Consider a Upper Triangular matrix U of dimension n, U x = b $$\begin{bmatrix} u_{1,1}&u_{1,2}&\dots&u_{1,n}\\ 0&u_{2,2}&\dots&u_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&u_{n,n} \end{bmatrix} \begin{bmatrix} x_{1}\\ x_{2}\\ \vdots\\ x_{n} \end{bmatrix} = \begin{bmatrix} b_{1}\\ b_{2}\\ \vdots\\ b_{n} \end{bmatrix}$$ Solving Procedure: $$x_{n} = \frac{b_{n}}{u_{n,n}}$$ Moving to ${n-1}^{th}$ row, there are two coefficients related to $x_{n}$(already known) and $x_{n-1}$ (unkown). $$u_{n-1,n-1}x_{n-1}+u_{n-1,n}x_{n}=b_{n-1}$$ Solving for unknown $x_{n}$, we have $$x_{n-1} = \frac{b_{n-1}-u_{n-1,n}x_{n}}{u_{n-1,n-1}}$$ Generalizing for $i^{th}$ row, $$x_{i} = \frac{1}{u_{i,i}}\left( b_{i}-\sum_{j=i+1}^{n}u_{i,j}x_{j} \right)$$ A ...
Read more

Solve Poissons Equation using 5-Point Stencil with SOR method

Image
Solution Poisson's Equation In a Rectangular Domain Using 5-Point Stencil Poisson's equation in a square domain $0 \leq x,y \leq 1$ with source term given by $q(x)=100*\sin(\pi x)*\sin(\pi y)$. The boundaries of the domain are maintained at $u=0$. Determine the solution using finite difference method. Solution: Discretizing domain uniformly with step size h and k along x and y direction. Using five point stencil, equation at each node is given by $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}+\frac{u_{i,j+1}-2u{i,j}+u{i,j-1}}{k^{2}}=-q(i,j)$$ $$u_{i,j}=\frac{1}{\frac{2}{h^{2}}+\frac{2}{k^{2}}}(\frac{u_{i+1,j}+u_{i-1,j}}{h^{2}}+\frac{u_{i,j+1}+u{i,j-1}}{k^{2}}+q(i,j))$$ where $q(i,j)=100*\sin(\pi x_{i})*\sin(\pi y_{j})$ Here a MATLAB code is given below. Choosen $m=n=11$ grid to solve the problem using SOR (successive overrelaxation). Mathematica Code: In this algo...
Read more