I would like to use the approximation of $f=a_x|x|+a_y|y|$ in my finite-element procedure, where $a_x$ and $a_y$ are constants and $x$ and $y$ are the unknowns that are subject to minimization in the finite-element procedure. So the first and second differentiations are needed for computing the residual and the tangent matrices. For the first trial, I used $\tilde{f}=\sqrt{a_x x^2 + a_y y^2 +\rho^2}$ with a reasonably small $\rho$. However, it has many convergence difficulties when $\rho$ is small, i.e. the corresponding element in the tangent matrix drastically increases when $x$ or $y$ or both are very small. On the other hand, I don't want to implement big values for $\rho$, because it affects the results.
Has anyone had some similar experiences with this kind of issues? Are there any better functions that can approximate my $f$ function in a more sophisticated manner for finite-element implementation? I should add here that when my $f$ function is just a function of one variable $x$ then the square root is fine for me. However, the problem arises when I have two variables.
I appreciate any help or suggestion.
Thanks in advance.