HW #1
For N=5, the node map looks like the figure below. The nodes start from 0 at the bottom left corner.
Sketched solution map of analytic solution:
The solution map of U has some clearly defined contour lines at theta = pi/6 and theta = pi/2, because of the cos(3*theta) term. At theta = pi/3 and theta = 0, we have the max and min values of the solution.
Solution Process:
I next created arrays for x(k), y(k), r(k), and theta(k). They are saved here:
x.dat
y.dat
r.dat
theta.dat
After verifying the analytic solution, which I did by plugging in the solution U into Laplace's Equation with cylindrical coordinates, I then took center difference equations for each term.
Next I created 2D contours for the analytic solution to compare to the initial solution. For N=5, the contour looked like the following figure.
Then, continuing on to the numerical solution, I used Matlab's left division operator to directly solve Ax = B. The A and B matrices are shown below for N=5.
For N=5, the solution U and the current vector, grad(U).
For N=10, the solution and current vector.
It's hard to see from the N=5 solutions whether or not the numerical solution agrees with the initial solution map of U. Since the mesh with 5 nodes, including the boundaries, is so coarse, it's difficult to accurately assess the PDE across the geometry. But as N is increased to 10, the general picture of the solution becomes clearer, and the similarity to the solution map can be seen.
Absolute Error:
Once I had both NxN matrices for the analytic solution and the numerical solution, it was a simple step to calculate the RMS error. The error map is shown on the geometry for both N=5 and N=10.
They both show that the error is very close to 0 for the majority of the geometry, but as N increases, there is a corner where the error begins to look problematic. This could be do to an issue with the molecule along that boundary.
RMS Error vs ∆r:
In my simulation, I used the single precision option in MATLAB, which still contained more than six significant figures. I found that as N increased, the error decreased as expected, and went as O(h^2).
Another expected result was that the error continued to decrease as N increased, but hit a saturation point and started to increase thereafter. This would be the point where the round off error begins to dominate the error calculation.
“Converged” (high N) :
For high N, I chose N=125 to show the converged numerical solution for the analytic solution and the current vector. With the mesh, it becomes a little difficult to see the solution, but the colormap agrees with the initial sketch.
For "high" N, this would be at the point that where round-off starts to dominant the error calculation. If my solution had worked correctly, I would have seen an uptick in the error, which would signal N has reached a sufficient level.