When using float16, computations are carried out in float16 precision. I would like to apply this to an explicit finite difference scheme for the heat equation, using

$$\mathrm{fl}(x \pm y) = x(1+\alpha) \pm y(1+\beta), \qquad |\alpha|, |\beta| \leq u$$

$$\mathrm{fl}(x \mathbin{\mathrm{op}} y) = (x \mathbin{\mathrm{op}} y)(1+\delta), \qquad |\delta| \leq u, \qquad \mathrm{op} \in {*, /}$$

I need to be sure which method is correct. The more rigorous method is the second one, because it follows the actual order of floating-point operations.

The first expression,

$$\tilde{u}_i^{,n+1} = \tilde{u}i^{,n} + r\Big(\tilde{u}{i+1}^{,n}(1+\varepsilon_1) - 2\tilde{u}i^{,n}(1+\varepsilon_2) + \tilde{u}{i-1}^{,n}(1+\varepsilon_3)\Big)(1+\varepsilon_4)$$

is a global model. It is useful for a qualitative analysis of round-off errors, but it does not describe exactly how the computer performs the computation.

The second expression,

$$t_1 = \mathrm{fl}(2u_i^n), \qquad t_2 = \mathrm{fl}(u_{i+1}^n - t_1), \qquad t_3 = \mathrm{fl}(t_2 + u_{i-1}^n), \qquad t_4 = \mathrm{fl}(r,t_3), \qquad u_i^{n+1,\mathrm{fl}} = \mathrm{fl}(u_i^n + t_4)$$

is more accurate because it describes the propagation of round-off errors according to the actual sequence of operations executed by the code.

I need this in order to prove the stability of the scheme using the von Neumann method, but I am really lost. Since this is my first time dealing with this, I cannot find the stability condition, the numerical value of the CFL condition, or how round-off errors influence the result.