Inverse potential problem

The inverse potential problem seeks the shape of a homogeneous source from measurements on a circle enclosing it. The operator maps the shape of the heat source to the heat flux measured on that circle. The heat distribution satisfies

\[\begin{split}\begin{cases} \Delta u = 1_K & \text{ in } \Omega \\ u = 0 & \text{ on } \partial\Omega \end{cases}\end{split}\]

where \(\partial\Omega\) is the measurement circle and \(K\) is the heat source. The forward operator maps the shape of the heat source to the Neumann data:

\[\partial K \mapsto \frac{\partial u}{\partial\nu}|_{\partial\Omega}.\]

References:

• 6. Hettlich & W. Rundell “Iterative methods for the reconstruction of an inverse potential problem”, Inverse Problems, 12 (1996) 251–266.

• 20. Hohage “Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and an inverse scattering problem”, Inverse Problems, 13 (1997) 1279–1299.

Import auxiliary libraries and define logging

We use numpy arrays as the main data structure for the domain and rely on matplotlib.pyplot as the tool for plots.

To observe the information of the solver while computing we define the logging string and set the log level to info.

import logging

import matplotlib.pyplot as plt
import numpy as np
from potential import Potential


logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s %(levelname)s %(name)-40s :: %(message)s'
)


example_dir = "../../../../examples/potential/"
Path(example_dir).resolve()
sys.path.insert(0, str(example_dir))

Define Operator and regularization setting

The regularization setting uses a Sobolev penalty and an \(L^2\) data-fidelity term. Here the abstract Sobolev space defaults to smoothness index 1.

from regpy.solvers import Setting
from regpy.hilbert import L2, Sobolev


#Forward operator
op = Potential(
    radius=1.3
)

setting = Setting(op=op, penalty=Sobolev, data_fid=L2)

2026-01-29 15:09:39,191 WARNING  Setting              :: Setting does not contain any explicit data.
2026-01-29 15:09:39,191 WARNING Setting                                  :: Setting does not contain any explicit data.

The exact solution and data

For the exact solution, we define a star-shaped domain with boundary $ \sqrt{3\cos(t)^2+1}/2$. We then compute exact data and add noise with a relative level of 3%.

#Exact data and Poission data
exact_solution = op.domain.radialfct2curve(lambda t: np.sqrt(6*np.cos(1.5*t)**2+1)/3)
exact_data = op(exact_solution.radial())
noise = op.codomain.randn()
noise = 0.03*setting.h_codomain.norm(exact_data)/setting.h_codomain.norm(noise) * noise
data = exact_data + noise

fig = plt.figure(figsize=(7,7))
plt.plot(*exact_solution.z)
ax = plt.gca()
ax.set_xlim([-1.1, 1.1])
ax.set_ylim([-1.1, 1.1])
plt.show()

Solving the problem

With an initial guess that is the circle of radius 1 we use a Newton CG solver to get the solution. For a reference we plot the initial guess with the true objective.

We combine the discrepancy principle with a maximum of 100 iterations to limit the run time.

#Initial guess
init = op.domain.circle(radius=0.2)

fig = plt.figure(figsize=(7,7))
plt.plot(*exact_solution.z)
plt.plot(*init.z)
ax = plt.gca()
ax.set_xlim([-1.1, 1.1])
ax.set_ylim([-1.1, 1.1])
plt.show()

from regpy.stoprules import CountIterations, Discrepancy

from regpy.solvers.nonlinear import NewtonCG
solver = NewtonCG(
    setting,
    data,
    init = init.coeff,
    cgmaxit=50,
    rho=0.6
)

stoprule = (
    CountIterations(100) +
    Discrepancy(
        setting.h_codomain.norm, data,
        noiselevel = setting.h_codomain.norm(noise),
        tau=1.1
    )
)

2026-01-29 15:09:39,825 WARNING  Discrepancy          :: Initialization of Discrepancy with two positional arguments (norm, data) deprecated. Use one positional argument (noiselevel) and specifiy norm and data via a keyword argument setting!
2026-01-29 15:09:39,825 WARNING Discrepancy                              :: Initialization of Discrepancy with two positional arguments (norm, data) deprecated. Use one positional argument (noiselevel) and specifiy norm and data via a keyword argument setting!

Run the iterative regularization method and plot results in each step

for reco,reco_data in solver.until(stoprule=stoprule):
    fig, axs = plt.subplots(1, 2,figsize=(15,7))
    axs[0].set_title('Obstacle')
    axs[0].set_xlim([-1.1, 1.1])
    axs[0].set_ylim([-1.1, 1.1])
    axs[1].set_title('Heat flux')

    axs[0].plot(*exact_solution.z)
    axs[0].plot(*op.curve.z)

    axs[1].plot(exact_data.real, label='exact')
    axs[1].plot(reco_data.real, label='reco')
    axs[1].plot(data.real, label='measured')
    axs[1].legend()
    axs[1].set_ylim(ymin=0)

    plt.show()

2026-01-29 15:09:40,045 INFO     CombineRules         :: it. 0>=100 | discr./noiselevel = 77.03< 1.10
2026-01-29 15:09:40,045 INFO CombineRules                             :: it. 0>=100 | discr./noiselevel = 77.03< 1.10

2026-01-29 15:09:40,292 INFO     CombineRules         :: it. 1>=100 | discr./noiselevel = 13.67< 1.10
2026-01-29 15:09:40,292 INFO CombineRules                             :: it. 1>=100 | discr./noiselevel = 13.67< 1.10

2026-01-29 15:09:40,490 INFO     CombineRules         :: it. 2>=100 | discr./noiselevel = 2.74< 1.10
2026-01-29 15:09:40,490 INFO CombineRules                             :: it. 2>=100 | discr./noiselevel = 2.74< 1.10

2026-01-29 15:09:40,684 INFO     CombineRules         :: it. 3>=100 | discr./noiselevel = 1.81< 1.10
2026-01-29 15:09:40,684 INFO CombineRules                             :: it. 3>=100 | discr./noiselevel = 1.81< 1.10

2026-01-29 15:09:40,958 INFO     CombineRules         :: it. 4>=100 | discr./noiselevel = 1.28< 1.10
2026-01-29 15:09:40,958 INFO CombineRules                             :: it. 4>=100 | discr./noiselevel = 1.28< 1.10

2026-01-29 15:09:41,230 INFO     CombineRules         :: it. 5>=100 | discr./noiselevel = 1.00< 1.10
Rule Discrepancy(noiselevel=0.012891189621064078, tau=1.1) triggered.
2026-01-29 15:09:41,230 INFO CombineRules                             :: it. 5>=100 | discr./noiselevel = 1.00< 1.10
Rule Discrepancy(noiselevel=0.012891189621064078, tau=1.1) triggered.