Implementing Edge Preserving Diffusion (Anisotropic Diffusion)

Question

I am trying to implement edge preserving diffusion. Recall the general diffusion equation to be:

$$\DeclareMathOperator{\Div}{div}\delta_t u = \Div(g \nabla u)$$

Where $g$ is the speed of the diffusion process and $u$ is the image. Now if we want an edge preserving diffusion, we can use a non linear diffusion speed as in Scale Space and Edge Detection Using Anisotropic Diffusion:

$$g(|\nabla u|) = \frac{1}{\sqrt{1+\frac{|\nabla u|}{\lambda^2}}}$$

To implement this smoothing effect, we can simply rewrite the diffusion equation as:

$$\delta_t u = \delta_x(g(|\nabla u|)\delta_xu) + \delta_y(g(|\nabla u|)\delta_yu)$$

And then perform gradient descent: $$u_{t+1} = u_t + \tau\delta_t u$$

This is implemented as follows:

lam = 0.1
im = smisc.imread('lena.png', mode='L')
for i in range(100):
    ux = cv2.Sobel(im, cv2.CV_64F, 1, 0, ksize=5)  # x
    uy = cv2.Sobel(im, cv2.CV_64F, 0, 1, ksize=5)  # y

    g = 1 / (1 + (ux ** 2 + uy ** 2) / lam ** 2) ** 0.5
    utx = cv2.Sobel(g * ux, cv2.CV_64F, 1, 0, ksize=5)  # x
    uty = cv2.Sobel(g * uy, cv2.CV_64F, 0, 1, ksize=5)  # y

    im = im + (utx + uty) * 0.01  # gradient descent
plt.imshow(im, cmap='gray', vmin=0, vmax=255)
plt.show()

This works well and smooths the image while preserving the edges. I decided that I would also try to implement this a different way, in particular I would solve the equations until the end and implement the smoothing. In particular I did: $$\begin{align} \delta_t u &= \Div(g \nabla u) \\ &= \delta_x(g(|\nabla u|)\delta_xu) + \delta_y(g(|\nabla u|)\delta_yu) \\ &= u_{xx}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{1}{2}} - \frac{1}{\lambda^2}u_x^2u_{xx}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{3}{2}} \\ &+ u_{yy}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{1}{2}} - \frac{1}{\lambda^2}u_y^2u_{yy}\big(1 + \frac{u_x^2 + u_y^2}{\lambda^2}\big)^{-\frac{3}{2}} \\ \end{align}$$

Using gradient descent this can be implemented as:

lam = 0.1
im = smisc.imread('lena.png', mode='L')
for i in range(100):
    ux = cv2.Sobel(im, cv2.CV_64F, 1, 0, ksize=5)  # x
    uxx = cv2.Sobel(ux, cv2.CV_64F, 1, 0, ksize=5)  # x
    uy = cv2.Sobel(im, cv2.CV_64F, 0, 1, ksize=5)  # y
    uyy = cv2.Sobel(uy, cv2.CV_64F, 0, 1, ksize=5)  # y

    utx = uxx * (1 + (ux**2 + uy**2) / lam**2) ** -0.5 - 1 / lam**2 * ux**2 * uxx * (1 + (ux**2 + uy**2) / lam**2) ** (-1.5)
    uty = uyy * (1 + (ux**2 + uy**2) / lam**2) ** -0.5 - 1 / lam**2 * uy**2 * uyy * (1 + (ux**2 + uy**2) / lam**2) ** (-1.5)

    im = im + (utx + uty) * 0.01 # gradient descent
plt.imshow(im, cmap='gray', vmin=0, vmax=255)
plt.show()

Interestingly this smooths as well, however, there appears to be some weird dot artifacts appearing in some locations, does anyone know why this happens and what is wrong with the second implementations? I show the outputs below:

Original image:

Implementation 1:

Implementation 2: