I am attempting to reproduce published results of a finite differencing problem in space. The original author uses Matlab and I am using Python.
In my problem, I have a diagonal system of equations: $$ \text{diag}(A)x = b $$ Where $A\in \mathbb{R}^{N\times N}$ and $x$ and $b\in \mathbb{R}^{N\times 1}$.
Because A is diagonal, I have set it up to be simply a series of division operations: $$ x_1 = \frac{b_1}{A_{1,1}} \\ \vdots \\ x_N = \frac{b_N}{A_{N,N}} $$
So far the problem is rather straightforward. However, $b_i$ and $A_{i,i}$ happen to be very small values - on the order of $10^{-13}$ and $10^{-19}$, respectively, thus forcing $x_i$ to be on the order of $10^6$.
Because of the magnitudes of the input of the problem, the 2-norm of my calculated data with the published data is rather interesting. For example, for $i=1$, the difference between my data and the authors for $b_1$, $A_{1,1}$, and $x_1$ is:
$1.677\times 10^{-22}$, $1.3979\times 1.379210^{-28}$, and $8.092\times 10^{-03}$, respectively.
I am rather surprised how the differences between the inputs can be so low (within machine precision) yet the output error is so relatively high.
Can someone explain why this is occurring?
We can rewrite the equation as \begin{equation} \ln(x) = \ln(b)-\ln(A) \end{equation} Then suppose we compare two roots $x_1$ and $x_2$, then we can have \begin{equation} \ln(x_1)-\ln(x_2) = \ln(b_1)-\ln(A_1)-(\ln(b_2)-\ln(A_2)) \end{equation} This leads to \begin{equation} \ln(x_1/x_2) = \ln(b_1/b_2)+\ln(A_2/A_1) \end{equation} As you mentioned that \begin{equation} b_1 = b_2+1.677\times 10^{-22} \end{equation} \begin{equation} A_2 = A_1+1.3979\times 10^{-28} \end{equation} Then we can have that \begin{equation} \ln(x_1/x_2) = \ln(1+1.677\times 10^{-22}) + \ln(1+1.3979\times 10^{-28}) \end{equation} We also know that if $x<1$ then we can have \begin{equation} \ln(1+x) = x-1/2x^2 +1/3x^3+ .... \end{equation} Considering the machine eps is quite small, so we can know that \begin{equation} \ln(x_1/x_2) \approx 1.677\times 10^{-9} + 1.3797\times 10^{-9} \end{equation} So we can have \begin{equation} x_1/x_2 \approx e^{10^{-9}} \end{equation} We also know that if $x<1$ we can have \begin{equation} e^{x} = 1+x + 1/2 x^2 + 1/6x^3 +... \end{equation} So we can have that \begin{equation} x_1/x_2 \approx 1+ 10^{-9} \end{equation} Then we can compare $x_1$ and $x_2$ as \begin{equation} x_1 -x_2 = (x_1/x_2-1) \times x_2 \approx 10^{-9}*x2 \end{equation} As you mentioned that $x_2$ is around $10^6$, so we can know that $x_1-x_2\approx 10^{-3}$
It is worth noting that this is just a rough calculation, but it will give you a sense for why the difference could be much larger.