1
$\begingroup$

In my calculation (of a simple heat equation, for testing) using the Newton method, I tried to replace the full Jacobian matrix with an approximation vector, i.e. replacing $J$ in $$J(u)\delta u=-F(u)$$ with $$J\vec{v}\approx\frac{F(\vec{u}+\varepsilon \vec{v}) - F(\vec{u})}{\varepsilon}$$ which should be possible, after I solve the system using a GMRES-solver. Here $\varepsilon$ is calculated using $$\varepsilon=\frac{1}{n\vert\vert\vec{v}\vert\vert_2}\sum_{i=1}^n\left(b\vert u_i\vert+b\right)\text{, }b=10^{-6}$$ $F(u)$ is the function itself, i.e. $$F(u)=\nabla^2u+f$$ Now I encountered two problems:

  • Usually in the first iteration $\vert\vert\vec{v}\vert\vert_2=0$, resulting in infinity for $\varepsilon$. My current solution is to set $\varepsilon=10^{-6}$ in that case, but is that correct? I could not find any solution for that in the literature
    • When considering the value for $F(u)$ at each step, it decreases (as it should) when using the full version of $J$, but increases when using the approximation. Is there a way how can I narrow down the possible problems here?

How can I solve those problems (or is there literature about them which I could not find yet)?

Improve this question
$\endgroup$
3
  • $\begingroup$ GMRES does not need the matrix itself, only a response in the form of a matrix-vector product. This matrix-vector product can be approximated using a finite difference resulting in a vector (also called the jacobian-free newton-krylov method). $\endgroup$
    – arc_lupus
    Dec 12, 2018 at 9:53
  • 2
    $\begingroup$ You have a linear PDE. Why would you compute the application of the matrix via finite differences when you know what the exact value is? $\endgroup$
    – Wolfgang Bangerth
    Dec 12, 2018 at 21:57
  • $\begingroup$ @WolfgangBangerth: Because I use that equation as test case for the solver. My final application contains a non-linear system of equations. Before applying the method to the more complex stuff, I would like to calculate something with a known solution. $\endgroup$
    – arc_lupus
    Dec 13, 2018 at 7:08

2 Answers 2

Reset to default
6
$\begingroup$

Not sure where you get your equation for $\epsilon$, but ultimately your approximation for the Jacobian matvec operation is a finite difference approximation to the directed derivative of $F(\cdot)$. This means that you will want as accurate of an approximation to this directed derivative as possible while being sure to avoid numerical round-off issues.

For this situation, I have personally used the following method referred to as the Complex-Step Derivative Approximation:

\begin{align} \nabla F(\boldsymbol{x}) \cdot \boldsymbol{v} = \frac{\text{Im}(F(\boldsymbol{x} + i \epsilon \boldsymbol{v}))}{\epsilon} + O(\epsilon^2) \end{align}

where $\text{Im}(\cdot)$ maps a complex vector in $\mathbb{C}^n$ to $\mathbb{R}^n$ by just taking the imaginary part of each component. Note that if $\boldsymbol{v}$ is not a unit vector, you should compute the above quantity with the direction of $\boldsymbol{v}$ and then multiply the result by $\left\lVert \boldsymbol{v}\right\rVert$. The obvious thing to note here is that the above assumes you can modify your code to allow for feeding in complex numbers into $F(\cdot)$. If you can do this, this approximation has great numeric precision since it does not have round-off issues related to subtractions, like typical Finite Difference approximations. You can readily choose $\epsilon = 10^{-8}$ and get close to machine precision.

Improve this answer
$\endgroup$
11
  • $\begingroup$ Unfortunately my code does not allow the usage of imaginary values directly. Of course I can implement that, but I would like to avoid doing that... $\endgroup$
    – arc_lupus
    Dec 12, 2018 at 16:02
  • 3
    $\begingroup$ @spektr It's called "Complex-Step Derivative Approximation". You might want to name the method if you are introducing advanced techniques. $\endgroup$
    – P. G.
    Dec 12, 2018 at 16:31
  • 1
    $\begingroup$ @P.G. I didn’t know the name, I read the paper about it years ago. I’ll update my response. $\endgroup$
    – spektr
    Dec 12, 2018 at 20:52
  • 1
    $\begingroup$ +1$\,$ I'm a huge advocate of the Complex Step Approximation, but there are some subtle caveats. Since you will most likely be using a language with built-in support for complex arithmetic, you must be careful about using its built-in functions. For example, in Julia X' is simply the transpose when working with real variables, but automagically becomes the hermitian conjugate when X is promoted to a complex variable. This will ruin the CSA approximation. Analytic functions such as exp(X), sin(X), X^2, etc are fine, but another common built-in to avoid is norm(X) $\endgroup$
    – greg
    Feb 15, 2022 at 15:28
  • 1
    $\begingroup$ However, norm(X) can be replaced by something like sqrt(sum(X.*X)) and X' by permutedims(X) $\endgroup$
    – greg
    Feb 15, 2022 at 15:28
1
$\begingroup$

I finally found the solution to that problem: In order to be able to reuse as much code as possible, I had one single function for calculating the right hand side of my problem, i.e. $-F(u)$, and for calculating the residual value. This function returned $-F(u)$ when calling, but I forgot to remove the minus-sign when using the same function for the approximation, thus resulting in $$J\vec{v}\approx\frac{-F(\vec{u}+\varepsilon \vec{v}) + F(\vec{u})}{\varepsilon}$$ After multiplying everything with $-1$, the approach works as expected.

Improve this answer
$\endgroup$

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.