Is the behaviour of the following filter plausible?

I've designed a filter based on this equation

$$E(u) = \frac{1}{2}\int_{\Omega} \left[(u-f)^2 + g(\lVert \nabla u \rVert)\right]dxdy$$

Euler-Lagrange equations in this case are (should be) given by

$$u - f - \frac{1}{2}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}g'(\lVert \nabla u\rVert)\right) = 0$$

For given $\lambda > 0$ if I chose $g(x) = e^{-x/\lambda}$ the EL equations are $$u - f + \frac{1}{2\lambda}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}e^{-\lVert \nabla u\rVert/\lambda}\right) = 0$$

Continuous gradient descent in this case is given by

$$u_t = f - u - \frac{1}{2\lambda}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}e^{-\lVert \nabla u\rVert/\lambda}\right)$$

Basically my idea would be... if the gradient magnitude is very high then $u = f$ (i.e. I'd like to preserve strong edges) otherwise (small variations) $u$ should be tending to some constant, or anyway smaller variations.

Now... This is my implementation in opencv

void g(InputArray input, OutputArray output) {
//cv::exp(-im,out);
float lambda = 0.0000001;
Mat rval, arg = input.getMat();
//cv::exp(-arg / lambda2, tmp);
cv::exp(-arg/lambda,rval);
output.assign(rval);
}

void dg(InputArray input, OutputArray output) {
Mat rval, tmp, arg = input.getMat();
g(arg,tmp);
float lambda = 0.0000001;
rval = (-1.0/lambda)*tmp;
output.assign(rval);
}

void evolveImage(const cv::Mat & input, cv::Mat & output, void (*g)(cv::InputArray, cv::OutputArray), void(*dg)(cv::InputArray, cv::OutputArray), const int & k_max) {
int X = input.cols;
int Y = input.rows;

Mat imfloat;
input.convertTo(imfloat, CV_32FC1);
Mat u;
imfloat.copyTo(u);
cv::Mat & f = imfloat;

double minval, maxval;

minMaxLoc(u, &minval, &maxval);

Mat dx, dy;
Mat gprime, a, b;
for (auto k = 0; k < k_max; ++k) {

Sobel(u, dx, CV_32F, 1, 0, 3); //Estimation differences in x
Sobel(u, dy, CV_32F, 0, 1, 3); //Estimation differences in y
minMaxLoc(dx, &minval, &maxval);
minMaxLoc(dy, &minval, &maxval);
a = gprime.mul(dx);
b = gprime.mul(dy);
Sobel(a, dx, CV_32F, 1, 0, 3);
Sobel(b, dy, CV_32F, 0, 1, 3);

u = f + 0.5f*(dx + dy);

}

minMaxLoc(u, &minval, &maxval);
u.convertTo(output, CV_8UC1, 255 / (maxval - minval), -minval * 255 / (maxval - minval));
}


I have this messy function that can be called by a main if you do cut and paste

int continuousGradientDescentExample(int argc, char **argv) {
if (argc != 2) {
std::cout << "Need filename..." << std::endl;
return 1;
}

namedWindow("Test", WINDOW_AUTOSIZE);
imshow("Test", im);
waitKey(0);

Mat imfloat;
im.convertTo(imfloat, CV_32FC1);
Mat u(imfloat.size(), imfloat.type());

int & X = im.cols;
int & Y = im.rows;

evolveImage(imfloat, u, g, dg, 50);

std::cout << type2str(u.type()) << std::endl;

Mat output(Size(2 * X, Y), CV_8UC1);
im.copyTo(output(Rect(0, 0, X, Y)));
u.convertTo(u, CV_8UC1);
u.copyTo(output(Rect(X, 0, X, Y)));

namedWindow("Original vs Anisotropic diffusion", WINDOW_AUTOSIZE);
imshow("Original vs Anisotropic diffusion", output);

waitKey(0);

return 1;
}


The code is working (i.e. doesn't crash, you can try it if you want). But the thing is no matter how I choose $\lambda$ in the evolveImage the matrices $a,b$ are always $0$ which implies the update of $u$ doesn't produce any effect. To me the equation makes sense, but I can't understand why there's no difference at all.

Can anyone give me an insight?