Inverting a diagonal (almost) linear relationship in a numerically stable way I have a relationship of the form
$$y=Dx$$
where $D$ is a diagonal, square matrix with  non-negative entries. Some of the diagonal elements of $D$ may be zero. However, since $D$ is calculated by a computer, I don't know which ones are "really" zero, and which ones are just very small. That is, in practice none of them are exactly zero, just -1.34534634e-16 or something, but I believe that in the "true" matrix, were I to calculate it by hand, that entry would be a zero. Since I know from the theory that the matrix has non-negative entries, I can assume that tiny negative values are "really" zero, but it's less clear with a small positive value like 1.0234234e-05 or something.
I need to calculate $\tilde x$, where:

*

*If $D_{i,i}>0$, then $\tilde x_i=x_i=D_{i,i}^{-1} y_i$.

*If $D_{i,i}=0$, then $\tilde x_i=a_i$, where $a$ is some otherwise unrelated, constant vector.

How can I do this? Notice that if $D_{i,i}$ is very small, the two rules will give wildly different results depending on whether I assume it's actually zero or not.
For context: I have a linear ODE $u'=Du$, and I want to compute $u$ from $u'$. The vector $a$ mentioned above is in fact $u(0)$: when an entry of $D$ is zero I know that the corresponding entry in $u(t)$ doesn't change over time, so I can obtain it from $u(0)$. The reason I want to compute $u$ from $u'$ is because I have another ODE $v'=Au$, so by rewriting this as something like $v'=AD^{-1}u'$ I can integrate it to compute $v(t)$. $D$ is the diagonalization of a symmetric semi-postive-definite real matrix, as computed by some linear algebra package.
 A: Let's integrate $v' = AD^{-1} u'$. We get:
$$
v(t) = v(0) + AD^{-1} \left[u(t) - u(0)\right]. \tag{*}
$$
Note that for $D_{ii}$ that are close to zero the diagonal entries $D_{ii}^{-1}$ (which are large) are multiplied by $u_i(t) - u_i(0)$ which are small. But this large-small approach does not give an answer what is the value of the limit
$$
\lim_{D_{ii} \to 0} \frac{u_i(t) - u_i(0)}{D_{ii}} = {?} 
$$
Let's derive the explicit expression for $v$ in a slightly other way:
$$
u' = Du \implies u = e^{tD} u(0)
$$
$$
v' = Au = Ae^{tD} u(0) \implies v = v(0) + A \left[\int_0^t e^{sD} ds\right] u(0)
$$
$$
\int_0^t e^{sD} ds = \int_0^t \sum_{k=0}^\infty \frac{s^k D^k}{k!} ds = 
\sum_{k=0}^\infty \frac{t^{k+1} D^k}{(k+1)!}.
$$
when $D$ is not singular the series converge to $D^{-1} (e^{tD} - I)$ which perfectly matches (*).
When $D$ is singular we can pull the first term out obtaining
$$
\int_0^t e^{sD} ds = t I + t\sum_{k=1}^\infty \frac{t^k D^k}{(k+1)!}
$$
For $D_{ii}$ that are close to zero, the sum would have simply $t$ in the $i$-th diagonal element.
So the final expression for $v(t)$ takes the following forms:
$$
v(t) = v(0) +  tAu(0) + tA \left[\sum_{k=1}^\infty \frac{t^k D^k}{(k+1)!}\right] u(0) = 
v(0) +  tA \left[\sum_{k=0}^\infty \frac{t^k D^k}{(k+1)!}\right] u(0).
$$
This gives you an answer for zero-diagonal entries, but requires summation for the rest.
The term $\sum_{k=0}^\infty \frac{t^k D^k}{(k+1)!}$ can be computed as matrix function $\phi(x) = \sum_{k=1}^\infty \frac{x^k}{(k+1)!} = \frac{e^x - 1}{x}$ applied to the matrix $tD$. Matlab has built in function funm for that based on this paper (beware of the paywall). Scipy also has scipy.linalg.funm.
A: Hint:
As to discriminate between actual zero and non-zero values of D, a practical way allowed on modern computers is to calculate D at two different levels of digits.  At increasing the digits, the actual zeros shall move towards 0, while small values shall remain relatively stable.
In reply to your comment, in my old Mupad CAS I can fix the number of floating point digits for the whole session (for whichever calculation, by putting DIGITS:=20 or else).
Now, if the matrix you are diagonalizing is input precisely (in terms of rationals or irrationals that can be computed at the desired accuracy) and if it is not ill-conditioned,
then if you get a D value of , say, $1 \cdot 10^{-16}$ with 20 digits,  changing to e.g. 25 digits you should get approx. $1 \cdot 10^{-21}$ if that is an actual zero.
If it is an actual $1 \cdot 10^{-21}$ it should remain approximately so.
If instead the starting matrix is input as rounded floating point numbers, the above considerations are not applicable.
And if the matrix is ill conditioned, then none of the D values are reliable.
