3
$\begingroup$

I would like to solve a system of differential equations $u+\nabla(\nabla\cdot u)=f$ or in more detail

$a+\partial_t^2a+\partial_t\partial_xb+\partial_t\partial_yc=f$

$b+\partial_x^2b+\partial_x\partial_ta+\partial_x\partial_yc=g$

$c+\partial_y^2c+\partial_y\partial_tb+\partial_y\partial_xc=h$

with periodic boundary conditions on a regular grid by the method of Finite Differences. $a=a(t,x)$, $b=b(t,x)$ and $c=c(t,x)$ are my unknowns and $f=f(t,x)$, $g=g(t,x)$ and $h=h(t,x)$ are known.

I discretized my domain in time ($t$) and space ($x,y$). What I already have now is a big equation system. The matrix has full rank. For a grid with only 3 nodes (=27 nodes overall, =81 equations) in each direction it looks like this

enter image description here

When increasing the number of nodes the system gets to big very soon for my computers memory. That's why I actually want the system matrix to be symmetric in order to use Cholesky factorization (hopefully) or something else. But I am not sure if there exists a numbering of nodes which would make the system matrix symmetric.

In literature the commonly used example is the Laplace equation, where they just use a symmetric stencil in order to approximate $\Delta$. In contrast to that I have to use one dimensional finite difference for the approximation of e.g. $\partial_x^2$ in a 3 dimensional space. Furthermore I am using one-sided differences for the mixed derivatives at boundaries, which might destroy symmetric. However I am not sure if this really destroys symmetry.

Only because I found this thread I've got the hope that there is a numbering of my nodes, which would make my system matrix symmetric. But actually I cannot imagine that, because a non-symmetric matrix cannot become a symmetric one by only permuting rows and columns in the same way.

So is there a way to get a symmetric system matrix for my problem? If not, what would you do instead in order to solve the system?

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ Have you tried using a preconditioned iterative linear solver? Aside from possibly the source term, you could use Jacobian-Free Newton-Krylov methods. If your source term is ill-conditioned, you might need the Jacobian matrix for that term, but all of the differencing operations can be implemented in matrix-free fashion. $\endgroup$
    – Geoff Oxberry
    Sep 30, 2014 at 21:18
  • $\begingroup$ The operator appears to be self-adjoint, so it should be possible to represent it as a symmetric matrix. This would be straightforward in a finite element framework; for a finite difference method, it's less immediately clear, but probably possible. $\endgroup$
    – Jesse Chan
    Oct 1, 2014 at 1:25
  • $\begingroup$ Your intuition is right: permuting rows and columns (i.e., choosing a different enumeration of degrees of freedom) can not make a matrix symmetric if it wasn't before. You'll have to choose a discretization that retains the symmetry from the continuous problem. $\endgroup$
    – Wolfgang Bangerth
    Oct 1, 2014 at 10:56
  • $\begingroup$ @Geoff Oxberry Yes, I tried it. Now I know from small sample problems, that I have to precondition, in order to solve it via a Krylov method (gmres in matlab). But unfortunately the incomplete LU factorization would need years to finish for convenient sizes of the problem. @Wolfgang Bangerth Thanks for the hint. Maybe I should replace my one-sided finite differences by symmetric once. I will try. $\endgroup$
    – Rob
    Oct 2, 2014 at 9:13
  • $\begingroup$ @Rob Cholesky is roughly twice as efficient as LU; if incomplete LU would need years to finish, incomplete Cholesky would also need years to finish. You could try something like AMG; AMG does well on pure diffusion problems, so it might work well here. $\endgroup$
    – Geoff Oxberry
    Oct 2, 2014 at 9:24

2 Answers 2

Reset to default
2
$\begingroup$

The differential operator is "grad after div" and so (in the continuous version) is symmetric and negative semi-definite. To get a symmetric discretization your discretization of the gradient has to be the (negative) adjoint of your discrete divergence. With finite differences this says, for example, that forward differences for the divergence imply backward differences for the gradient. Boundary conditions need special care.

Improve this answer
$\endgroup$
1
$\begingroup$

Can you use timestepping instead of discretizing in time? You have a second order time derivative, but you can rewrite that as a first order system, i.e. $\frac{\partial^2 u}{\partial t^2} = f$ can be rewritten by adding another equation

$\frac{\partial w}{\partial t} = f$

$\frac{\partial u}{\partial t} - w = 0,$

and then use time-stepping methods for first order systems ODEs. This way, you'd end up with only a 2D discretization and not a 3D one (which saves a ton of memory).

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks for this idea; maybe worth a look. I am not very familiar with time stepping methods. It's an iterative method, right? Actually I don't want to use an iterative method again, because the whole thing would already run in an loop, which optimizes an other problem. $\endgroup$
    – Rob
    Oct 2, 2014 at 9:22
  • $\begingroup$ It does run in a loop, but the memory requirements go from $O(N^3)$ to $O(N^2)$ if you use timestepping, and things are generally much faster. Am I correct in thinking that you are doing a 3D simulation? $\endgroup$
    – Jesse Chan
    Oct 2, 2014 at 21:32
  • $\begingroup$ No, it's just a sub problem which pops up in a bigger optimization problem. It is not directly connected to a situation in reality. $\endgroup$
    – Rob
    Oct 6, 2014 at 8:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.