-1

I am currently programming a code to find the equilibrium function that satisfies the poisson equation in 2D. In order to do this I use finite difference methods and the discrete equation I want to satisfy is: $$\frac{T_{j, i+1} + T_{j, i-1} + T_{j+1, i} + T_{j-1, i} - 4T_{j, i}}{h^2} = q_{j,i},$$ where $T_{j,i}$ is a temperature on a rectangular grid. In order to find the resulting array $\vec{T}$ that satisfies the equation I write it as a matrix equation: $$\begin{pmatrix} -4 & 1 & & & \\ 1&-4&1&\\&\ddots&\ddots&\ddots& \\ & &\ddots &\ddots& 1\\ & & & 1 & -4 \end{pmatrix}\begin{pmatrix}T_1\\T_2\\\vdots \\T_{N-1} \\T_{N} \\ \end{pmatrix} = \begin{pmatrix}q_1\\q_2\\\vdots \\q_{N-1} \\q_{N} \\ \end{pmatrix}, $$ where the array $\vec{q}$ is a constant array.

I iteratively solve this equation using the Gauss-Seidel method, however if I just run it forever and look at the average of $T$ after every iteration, it never converges to zero, i.e. if I run the iteration forever the average $T \rightarrow \infty$. Step-size by which the average temperature increases between iterations becomes constant eventually, but never goes the zero, how do I establish convergence of this method then, if not by a threshold on the change of the average temperature?

Does this imply my matrix is not convergent? I am fairly sure the matrix I use should be convergent, as it is widely claimed by many credible sources.

Akerai
  • 101
  • It's not clear to me from your description why the average $T$ should converge to zero -- it should converge to a finite value, however. If it doesn't, then you have implemented the iteration incorrectly. – Wolfgang Bangerth Dec 17 '17 at 21:22
  • Hi Akerai, I see that no answers have been accepted. I would provide an answer myself, but I would basically be copying information from the following two SE questions (full disclosure, I asked and answered them) https://scicomp.stackexchange.com/questions/21612/poisson-equation-finite-difference-with-pure-neumann-boundary-conditions and https://scicomp.stackexchange.com/questions/28593/why-naively-chopped-finite-difference-matrix-works-for-different-ode-boundary-co/28713#28713. – Charles Mar 20 '18 at 04:26
  • The gist of what you need to answer this question falls on two important things. First, where is your data located? The cell-center? Or the cell-corner? This actually changes the form of the boundary conditions (BCs) (which I discuss in the second link) which, as @user3209427 pointed out, is in fact very important. Second, what are your intended BCs? Dirichlet? Neumann? This is critical. While Dirichlet BCs should not cause any issues, pure Neumann BCs require delicate care. Pure Neumann BCs require that the net heat input to the system is zero from a discrete operator point of view. – Charles Mar 20 '18 at 04:33

1 Answers1

0

If you're continuously injecting heat into the system through the q function I can see the temperature continuously growing depending on your boundary conditions, but that's a formulation problem that you need to consider.

Halting iterative methods of solving equations is done using the residual:

If we have the system as:

$ A.x = b $

We define the nth residual as

$ r_n = (b - A.x_n) $

Then you just have to determine a stopping criteria, usually a norm (L1, L2, Linf) of the absolute or relative residual, where the relative residual is:

$r_{r,n} = r_n/r_i$

i is typically taken as zero, i.e. the first residual, but depending on a variety of factors you might instead choose to use the fifth or tenth residual after the system has settled a little.

user3209427
  • 407
  • 3
  • 10
  • This answer is not useful. The discretization is of the Laplace equation, so there is no time-dependent effect in play here. Furthermore, stopping criteria are also not useful: if the iteration diverges, then the iteration computes things the wrong way. – Wolfgang Bangerth Dec 17 '17 at 21:21
  • 1
    I didn't day anything about time-dependent effects, iterative methods have a starting point and if the boundary conditions are not set correctly (consider a fully insulated boundary) then the source term will cause a continuous growth in energy as the problem is ill-posed which is what Akerai is seeing. The rest of my post is relevant/useful because it actually answers the asked question on how to establish if the method has converged. – user3209427 Dec 18 '17 at 02:36
  • I read "if you're continuously injecting heat" as suggesting that there is a time trajectory. But regardless, I can't see anything wrong with the matrix as shown -- it clearly uses Dirichlet boundary conditions, for which the problem is perfectly well-posed. – Wolfgang Bangerth Dec 18 '17 at 04:15
  • With only three values per row in the matrix it way well be a well-posed problem in 1D but not 2D, and even then it appears that the Dirichlet condition hasn't been included in the RHS q vector which was why I considered it a problem with Nuemann boundary conditions, along with the infinity growth in energy. – user3209427 Dec 18 '17 at 04:37
  • 1
    Well, if there is no addition to the right hand side, then that would simply be zero Dirichlet boundary values. It's true that the original question only shows 3 of the 5 entries per row. But if there were Neumann boundary values, there would need to be rows that look differently than the ones shown. You're right that the stated form of the matrix cannot be correct, one way or the other -- I just assumed that the form as shown just omits the additional 2 off-diagonals. – Wolfgang Bangerth Dec 18 '17 at 14:12