In general Poissons equation refers to a class of partial di

In general, Poisson\'s equation refers to a class of partial differential equations like: nabla^2 phi = f where phi is the physical field (e.g. electric potential) and f is provided as a source term (e.g. electric charge distribution). Now we consider an axisymmetric (with respect to z-axis) problem. This means you can express the nabla^2 operator in its cylindrical coordinate form and let all Phi-derivatives be zero. The problem is simplified to two-dimensional with only rho and z dependence. To numerically solve the problem, the first step is to define the simulation domain: 0 lessthanorequalto rho lessthanorequalto R and 0 lessthanorequalto z lessthanorequalto H, where R and H are constants to be specified. Then we will specify the source term f(rho, z). To start with, you can set it zero everywhere. This is the case of Laplace equation. Next, specify the boundary condition for phi. You need first to think about what kind of boundary conditions are appropriate at the boundary of the simulation domain, i.e. when rho=0, rho=R, z=0 and z=H. Then you need to specify the boundary condition of phi at the boundary of a region inside the simulation domain, i.e. when G(rho, z)=0. The form of G(rho, z) is completely open. The purpose of this interior region is to include the cases when there are different materials in the simulation domain. Finally, choose a spatial step-length to discretize the simulation domain and use Finite Difference method to numerically solve the differential equation. Write on how you implemented these steps in MATLAB and how you validated the numerical simulation.

Solution

Entire question is about how do you solve Poisson\'s Equation in MATLAB.

Syntax

u = poisolv(b,p,e,t,f)

Description

u = poisolv(b,p,e,t,f) solves Poisson\'s equation with Dirichlet boundary conditions on a regular rectangular grid. A combination of sine transforms and tridiagonal solutions is used for increased performance.

The boundary conditions b must specify Dirichlet conditions for all boundary points.
The mesh p, e, and t must be a regular rectangular grid. For details on the mesh data representation, see initmesh.f gives the right-hand side of Poisson\'s equation.Apart from roundoff errors, the result should be the same as u = assempde(b,p,e,t,1,0,f).

Description

u = poicalc(f,h1,h2,n1,n2) calculates the solution of Poisson\'s equation for the interior points of an evenly spaced rectangular grid. The columns of u contain the solutions corresponding to the columns of the right-hand side f. h1 and h2 are the spacings in the first and second direction, and n1 and n2 are the number of points.

The number of rows in f must be n1*n2. If n1 and n2 are not given, the square root of the number of rows of f is assumed. If h1 and h2 are not given, they are assumed to be equal.

The ordering of the rows in u and f is the canonical ordering of interior points, as returned by poiindex.

The solution is obtained by sine transforms in the first direction and tridiagonal matrix solution in the second direction. n1 should be 1 less than a power of 2 for best performance.

0001% Numerical approximation to Poisson’s equation over the square [a,b]x[a,b] with
0002 % Dirichlet boundary conditions. Uses a uniform mesh with (n+2)x(n+2) total
0003 % points (i.e, n x n interior grid points).
0004 % Input:
0005 % pfunc : the RHS of poisson equation (i.e. the Laplacian of u).
0006 % bfunc : the boundary function representing the Dirichlet B.C.
0007 % a,b : the interval defining the square
0008 % n : n+2 is the number of points in either direction of the mesh.
0009 % Ouput:
0010 % u : the numerical solution of Poisson equation at the mesh points.
0011 % x,y : the uniform mesh.
0012 %
0013 function [u,x,y] = fd2poisson(pfunc,bfunc,a,b,n)
0014
0015 h = (b-a)/(n+1); % Mesh spacing
0016
0017 [x,y] = meshgrid(a:h:b); % Uniform mesh, including boundary points.
0018
0019 % Compute u on the boundary from the Dirichlet boundary condition
0020 ub = zeros(n,n);
0021 idx = 2:n+1;
0022 idy = 2:n+1;
0023 % West and East boundaries need special attention
0024 ub(:,1) = feval(bfunc,x(idx,1),y(idy,1)); % West Boundary
0025 ub(:,n) = feval(bfunc,x(idx,n+2),y(idy,n+2)); % East Boundary
0026 % Now the North and South boundaries
0027 ub(1,1:n) = ub(1,1:n) + feval(bfunc,x(1,idx),y(1,idy));
0028 ub(n,1:n) = ub(n,1:n) + feval(bfunc,x(n+2,idx),y(n+2,idy));
0029 % Convert ub to a vector using column reordering
0030 ub = (1/h^2)*reshape(ub,n*n,1);
0031
0032 % Evaluate the RHS of Poisson’s equation at the interior points.
0033 f = feval(pfunc,x(idx,idy),y(idx,idy));
0034 % Convert f to a vector using column reordering
0035 f = reshape(f,n*n,1);
0036
0037 % Create the D2x and D2y matrices
0038 z = [-2;1;zeros(n-2,1)];
0039 D2x = 1/h^2*kron(toeplitz(z,z),eye(n));
0040 D2y = 1/h^2*kron(eye(n),toeplitz(z,z));
0041
0042 % Solve the system
0043 u = (D2x + D2y)\\(f-ub);
0044 % Convert u from a column vector to a matrix to make it easier to work with
0045 % for plotting.
0046 u = reshape(u,n,n);
0047
0048 % Append on to u the boundary values from the Dirichlet condition.
0049 u = [[feval(bfunc,x(1,1:n+2),y(1,1:n+2))];...
0050 [[feval(bfunc,x(2:n+1,1),y(2:n+1,1))] u ...
0051 [feval(bfunc,x(2:n+1,n+2),y(2:n+1,n+2))]];...
0052 [feval(bfunc,x(n+2,1:n+2),y(n+2,1:n+2))]];