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 :(