Solution of symmeric/non-symmetric linear system

Question

I would like to understand what happens in the following:

I have a really simple Poisson problem, in 1D, with $u_0 = u_N = 0$. I assembled the stiffness matrix and the right-hand side, and I applied the BCs, then forced $A$ to be symmetric.

I'm studying the stability of an iterative method for solving this linear system while increasing $N$. I put $1e^{-10}$ as tolerance and everything went ok until $N =35000$ where around $1.2e^{-10}$ the residual starts to oscillate.

As a test I tried a solving in Matlab (using A\b), and also there the residual did not go below $1e^{-10}$.

So I removed the symmetry from the stiffness matrix, and I retried on Matlab, and now the residual is in the order of $1e^{-11}$.

So, is it possible that a simple manipulation for making the matrix symmetric could cause my solution to be "worse"?

Answer 1

It could have changed the condition number adversely, yes, which would make it harder to solve. How are you forcing this matrix to be symmetric? Should be said that if you're solving the linear problem 10 orders, you're probably doing pretty well.

As whpowell pointed out in their comment matlab will use a different solver routine for symmetric vs. asymmetric. It will choose between one of two solvers for a symmetric matrix and from a different two solvers for an asymmetric matrix. So you're not only changing the matrix condition number, you're changing the solver as well.

Answer 2

Not really an answer, but too long for a comment.

In my experience, the best way to preserve symmetry and deal with boundary conditions is to use just internal points. "Penalization methods" (or big number approach, I know it has several names) is another possibility .

If your enumeration is $x_1,\ldots,x_{n+1}$, then work just on $x_2, \ldots, x_{n}$, as the values of $u$ at $x_1,x_{n+1}$ are given (in your case, just $0$). Ig they're not $0$, they will move on the rhs-vector

You reduce the d.o.f of the stiffness matrix by $1$, and preserve the symmetry. You'll obtain a solution vector, after the lin. system solution process, consisting of just the $n-1$ interal values. You can then add the known values at the end in the sol. vector if you want.