Solving Poisson Equation with FFT in MATLAB

[rustpanjguh]

Hi,

As part of my homework, I wrote a MatLab code to solve a Poisson equation
Uxx +Uyy = F(x,y)

with periodic boundary condition in the Y direction and Neumann boundary condition in the X direction. I used the finite difference method in the X direction and FFT in the Y direction to numerically solve for Uxx and Uyy. Now the teacher is asking to use our numerical solver and test it for interesting functions.

My code solves the equations with above conditions and with smooth U with order two accuracy, and now the teacher wants us to run our code for interesting cases, like non-periodic functions or non-smooth functions.
Here are my questions:

1. I need help coming up with a periodic (in y) non-smooth function so I can test my code. Can someone help me find such functions? (MatLab implementation of function would be great)

2. When I find such function and run my code using it, then what should I expect? Should I expect that a 1st order approximation, or total explosion of solution?

3. As I mentioned For a non periodic function I got a 1st order error rate. Is that what I should expect?

 

[jvicens]

Hi there,

It's great that you have successfully implemented a numerical solver for the Poisson equation with periodic and Neumann boundary conditions. To test your code for interesting cases, you can try using non-periodic functions such as a Gaussian or a step function. These functions have different smoothness and can provide a good test for your code.

To implement a Gaussian function in Matlab, you can use the "normpdf" function which generates a normal probability distribution. You can then use this function to create a 2D Gaussian function with different means and standard deviations. For example, the following code creates a 2D Gaussian with a mean of 0 and a standard deviation of 1:

x = linspace(-5,5,100);
y = linspace(-5,5,100);
[X,Y] = meshgrid(x,y);
Z = normpdf(X,0,1).*normpdf(Y,0,1);
surf(X,Y,Z)

You can also try using a step function, which can be implemented in Matlab using the "heaviside" function. This will create a step function with different step locations and heights. For example, the following code creates a step function with a step at x=0 and a height of 1:

x = linspace(-5,5,100);
y = zeros(size(x));
y(x>=0) = 1;
plot(x,y)

When running your code with these functions, you should expect to see a solution that is not smooth and has discontinuities. Depending on the function and the accuracy of your code, you may see a first-order or even a second-order error rate. It is important to also check the convergence of your solution as you refine your grid size to ensure that your code is producing accurate results.

I hope this helps and good luck with your code testing! Let me know if you have any further questions.