5
$\begingroup$

Let $n,m\in \mathbb{N}$ be such that $m\ge n$.

Let $M_1\in \mathbb{R}^{n\times n}$, $\{M_{2},M_{3} \} \subset \mathbb{R}^{m\times m}$ be symmetric positive definite and computationally cheap to invert. You can think of them as a diagonal matrix. Let $S\in \mathbb{R}^{n\times m}$ be full rank. Let $B\in \mathbb{R}^m$. I am interested in solving the following linear system over a range of $\omega \in \mathbb{R}$. $$ \begin{align*} \left\lbrack \begin{array}{cc} -i \omega M_1& S\\ -S^\mathrm{T} & -i\omega M_2+ M_3 \end{array} \right\rbrack \left\lbrack \begin{array}{c} U\\Q \end{array} \right\rbrack = \left\lbrack \begin{array}{c} 0\\B \end{array} \right\rbrack. \end{align*} $$

If $M_1$ and $M_2$ were both $0$ this would be a saddle point problem and if $M_3$ was zero this would be a Helmholtz problem. I am interested in solving this equation over a range of values of $\omega$, so reusing decompositions might be useful. This system comes from a Telegrapher's equation.

I am writing to inquire if there is a body of research on these types of problem. I know the research for the saddle point problems and Helmholtz problem is pretty extensive. Any ideas or suggestions are appreciated.

Edit.

This is what I am planning on doing, without getting a better idea from here.

Writing $U= U_R + i U_I$ and $Q= Q_R + i Q_I$ and rewriting the system by equating real and imaginary parts we get

$$ \left\lbrack \begin{array}{cccc} -\omega M_1 & 0 & 0 & S\\ 0 & \omega M_1& S&0\\ 0& -S^\mathrm{T}& -\omega M_2 & M_3\\ -S^{\mathrm{T}}&0&M_3& \omega M_2 \end{array} \right\rbrack \left\lbrack \begin{array}{c} U_R\\ U_I\\ Q_R\\ Q_I \end{array} \right\rbrack = \left\lbrack \begin{array}{c} 0\\ 0\\ 0\\ B \end{array} \right\rbrack. $$

Eliminating $U_I$ and $U_R$ from these we get

$$ \left\lbrack \begin{array}{cc} \omega M_3 & S^\mathrm{T} M_1^{-1} S - \omega^2 M_2\\ S^\mathrm{T} M_1^{-1} S - \omega^2 M_2& \omega M_3 \end{array} \right\rbrack \left\lbrack \begin{array}{c} Q_R\\ Q_I \end{array} \right\rbrack = \left\lbrack \begin{array}{c} \omega B\\0 \end{array} \right\rbrack. $$

Now we have a standard saddle point problem. I plan on solving it using a Schur complement (inner/outer CG solves) method. Note that I will use the $\omega M_3$ matrices as "pivots" because there exists $\omega$ for which $S^\mathrm{T} M_1^{-1} S - \omega^2 M_2$ is singular. Having obtained $Q_I$ and $Q_R$ we can use the eliminated equations to recover $U_I$ and $U_R$.

What I don't like about this approach is that there is no reuse. I mean the method doesn't achieve any economies of scale even though it will solve many closely related problems.

Improve this question
$\endgroup$
10
  • $\begingroup$ How big are $n$ and $m$ likely to be? Have you tried a block factorization? That should let you reuse factorizations of all of the submatrices, though you will have to refactor $i\omega M_2+M_3$ each time you change $\omega$ unless there's some very special structure in $M_2$ and $M_3$. $\endgroup$
    – Bill Barth
    Oct 1, 2015 at 21:14
  • $\begingroup$ The matrices will be as can big as can fit on the computer's RAM. (normal workstation not a super computer). The $M$ matrices are actually very sparse and almost diagonal so inverting them for a given RHS is a trivial operation. $m$ will be $3$ times $n$ for $3$- dimensional problems and and $2$ times $n$ for $2$ dimensional problems. As for as factorizations, I really only think it would be worthwhile to to factor $S$, maybe into some sort of sparse SVD. $\endgroup$
    – fred
    Oct 1, 2015 at 21:49
  • 2
    $\begingroup$ Do block-wise Gaussian elimination. $\endgroup$
    – Bill Barth
    Oct 2, 2015 at 2:52
  • 2
    $\begingroup$ "Do block-wise Gaussian elimination" -- which is also known as forming the "Schur complement" in this context. $\endgroup$
    – Wolfgang Bangerth
    Oct 2, 2015 at 2:53
  • 1
    $\begingroup$ @BillBarth So nothing glaringly obviously wrong or inappropriate as far as you see it? Thanks. $\endgroup$
    – fred
    Oct 2, 2015 at 14:33

1 Answer 1

Reset to default
1
$\begingroup$

Sorry this is late, but I've been thinking about an iterative solution of this which might be an alternative. So you're solving $(\omega A+B)x = b$ for $x$ with $A$ easily invertible. Can you make an iterative version as $$\omega A x^{(i+1)} + Bx^{(i)} = b,$$ or equivalently $$x^{(i+1)} = \omega^{-1}A^{-1}b - \omega^{-1}A^{-1}Bx^{(i)}.$$ You can precalculate $A^{-1}b$ and $A^{-1}B$, and as @BillBarth said, you could try the $x$ for the previous $\omega$ as your starting guess. Maybe I'm missing something.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Was going to comment this but not enough rep. Can delete if it's useless. $\endgroup$
    – Steve
    Oct 5, 2015 at 13:37
  • $\begingroup$ Hi, I think this is valid. I will check and try to get an idea of the spectrum of the $A^{-1}B$ operator to see how fast this will converge. It looks like for large $\omega$ it will converge fast. $\endgroup$
    – fred
    Oct 5, 2015 at 14:03
  • $\begingroup$ Also, since you're only doing the inversion once, could check out $$x^{(i+1)} = B^{-1}b - \omega B^{-1}Ax^{(i)}$$ $\endgroup$
    – Steve
    Oct 5, 2015 at 14:36
  • $\begingroup$ This will only work for $\omega$ large enough. For your first example, we need that $\left\|\omega^{-1} A^{-1}B \right\| <1$. This is a pretty drastic restriction on $\omega$ in my application since $ \left\| A^{-1}B \right\| \to \infty$ as the mesh gets smaller. There might be a way to fix this though. Thanks for the idea. $\endgroup$
    – fred
    Oct 5, 2015 at 15:04

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.