How can I solve a non-linear functional problem numerically?

[Irid]

Hello,
I'm not really sure where does this question fit and what title should it bear, but here is my problem:

$\psi(x) \exp (a\psi(x)^2) = C f(x)$

given a positive definite f(x), find ψ(x) and the constant C, subject to the condition

$\int \psi(x)\, dx = 1$

I want to solve this numerically but I'm having a hard time devising a suitable algorithm... is this even a well-posed problem?

 

[economicsnerd]

If $a\geq 0$, then it's well-posed.

Indeed, $ye^{ay^2}$ is a strictly increasing function of $y$ which vanishes at $y=0$ and blows up as $y\to\infty$, so for each $C\geq 0$, there's a unique choice of function $\psi_C$ that satisfies your functional equation. Also, the solution $\psi_C$ is going to be as smooth as $f$ is; in particular, it's integrable if $f$ is (say) piecewise continuous. Moreover, $\psi_C$ increases (smoothly enough to make $\int\psi_C$ continuous in $C$ if it's finite-valued) with $C$, and you can use a bounding argument to show that $\psi_C$ can be made arbitrarily big or small as you vary $C$. It follows that as long as $f$ is such that, for small enough $C$, the solution $\psi_C$ has bounded integral (which comes for free if your domain of integration is bounded), there's a unique choice of $\psi, C\geq 0$ that solves your equation.

As for how to solve numerically, I have no idea, but hopefully getting a sense of the structure of your equation can be helpful...

 

[Irid]

If $a\geq 0$, then it's well-posed.

Indeed, $ye^{ay^2}$ is a strictly increasing function of $y$ which vanishes at $y=0$ and blows up as $y\to\infty$, so for each $C\geq 0$, there's a unique choice of function $\psi_C$ that satisfies your functional equation. Also, the solution $\psi_C$ is going to be as smooth as $f$ is; in particular, it's integrable if $f$ is (say) piecewise continuous. Moreover, $\psi_C$ increases (smoothly enough to make $\int\psi_C$ continuous in $C$ if it's finite-valued) with $C$, and you can use a bounding argument to show that $\psi_C$ can be made arbitrarily big or small as you vary $C$. It follows that as long as $f$ is such that, for small enough $C$, the solution $\psi_C$ has bounded integral (which comes for free if your domain of integration is bounded), there's a unique choice of $\psi, C\geq 0$ that solves your equation.

As for how to solve numerically, I have no idea, but hopefully getting a sense of the structure of your equation can be helpful...

OK, great, that's some fancy math right there! Indeed, a>0 always. The f(x) will typically be known only numerically on the grid points and is not an elementary function, but generally it will resemble something between a Gaussian and an exponential decay, so it will have a maximum near the origin and will rapidly decay to zero at infinity, and certainly has a finite norm.

Actually, the problem is in 3D, so replace all instances of x by (x,y,z). My idea for the algorithm is as follows: take the log of the equation to obtain (the log will damp the errors)

$\log(\psi) + a\psi^2 = \log(C) + \log f$

Then guess some (normalized) function ψ1, say ψ1 = 1 (we are working in a finite box 1x1x1). Obviously this does not satisfy the equation, so we suppose that a better guess would be

$\psi_2 = \psi_1 + \Delta \psi_1$

where the added term is hopefully going to get small eventually. Then we Taylor expand the main equation to get

$\Delta\psi_1 = \frac{\psi_1}{1+2a\psi_1^2}(\log(Cf/\psi_1)-a\psi_1^2)$

Then repeat the last step until the delta gets as small as needed. I just don't see how to get the C out of this and how to ensure normalization...

 

[Irid]

Hmmm... If my initial guess is normalized, all I need to do is to ensure that the norm of Δψ is zero at each iteration. This can be naturally done by choosing a suitable log(C) = int f1/int f2.

I've just tried this, and the whole thing converges within 10 iterations max haha :D But apparently I need a very fine mesh, otherwise non-physical fluctuations show up :(