.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_ch3\ch3_4_3.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_ch3_ch3_4_3.py: ============================================= 3.4.3 Inverse Laplace operator ============================================= The inverse Laplace operator is a useful tool in many fluid mechanics, image processing, signal filtering etc. In the CodPy framework, this operator is implemented as the pseudo-inverse of the Laplacian matrix built from a kernel. Overview -------- The Laplacian operator $\Delta_k$ is defined in CodPy via kernel methods as a matrix $\Delta_k(X, Y) \in \mathbb{R}^{N_x \times N_x}$ constructed from a set of input points $X$ (and optionally $Y$). It generalizes the classical Laplacian to the setting of RKHS. The **inverse Laplacian** is simply the matrix inverse (or pseudo-inverse) of this operator: $$ \Delta_k^{-1}(X, Y) = \left( \Delta_k(X, Y) \right)^{-1} \in \mathbb{R}^{N_x \times N_x} $$ This operator provides a way to "undo" the effect of the Laplacian on a function. It is particularly useful when solving elliptic partial differential equations. CodPy Implementation -------------------- .. GENERATED FROM PYTHON SOURCE LINES 34-172 .. code-block:: Python # Importing necessary modules import os import sys from matplotlib import pyplot as plt curr_f = os.path.join(os.getcwd(), "codpy-book", "utils") sys.path.insert(0, curr_f) import numpy as np # from codpy.plotting import plot1D # Lets import multi_plot function from codpy utils from codpy.plot_utils import multi_plot # Define the sinusoidal function def periodic_fun(x): """ A sinusoidal function that generates a sum of sines based on the input ``x``. """ from math import pi sinss = np.cos(2 * x * pi) if x.ndim == 1: sinss = np.prod(sinss, axis=0) ress = np.sum(x, axis=0) else: sinss = np.prod(sinss, axis=1) ress = np.sum(x, axis=1) return ress + sinss def nabla_my_fun(x): from math import pi import numpy as np sinss = np.cos(2 * x * pi) if x.ndim == 1: sinss = np.prod(sinss, axis=0) D = len(x) out = np.ones((D)) def helper(d): out[d] += 2.0 * sinss * pi * np.sin(2 * x[d] * pi) / np.cos(2 * x[d] * pi) [helper(d) for d in range(0, D)] else: sinss = np.prod(sinss, axis=1) N = x.shape[0] D = x.shape[1] out = np.ones((N, D)) def helper(d): out[:, d] += ( 2.0 * sinss * pi * np.sin(2 * x[:, d] * pi) / np.cos(2 * x[:, d] * pi) ) [helper(d) for d in range(0, D)] return out # Function to generate periodic data def generate_periodic_data_cartesian(size_x, size_z, fun=None, nabla_fun=None): """ Generates 2D structured Cartesian grid data for x and z domains, and evaluates a given function and optionally its gradient. Parameters: - size_x: number of points per axis for x (grid will be size_x^2) - size_z: number of points per axis for z (grid will be size_z^2) - fun: function to evaluate at each point - nabla_fun: optional gradient function to evaluate Returns: - x, z: 2D Cartesian grids of shape (N, 2) - fx, fz: function values at x and z - nabla_fx, nabla_fz (if nabla_fun is provided) """ def cartesian_grid(size, box): lin = [np.linspace(box[0, d], box[1, d], size) for d in range(2)] X, Y = np.meshgrid(*lin) return np.stack([X.ravel(), Y.ravel()], axis=1) # Define domain boxes X_box = np.array([[-1, -1], [1, 1]]) Z_box = np.array([[-1.5, -1.5], [1.5, 1.5]]) # Generate Cartesian grids x = cartesian_grid(size_x, X_box) z = cartesian_grid(size_z, Z_box) # Function evaluations fx = fun(x).reshape(-1, 1) if fun else None fz = fun(z).reshape(-1, 1) if fun else None if nabla_fun: nabla_fx = nabla_fun(x) nabla_fz = nabla_fun(z) return x, fx, z, fz, nabla_fx, nabla_fz return x, fx, z, fz # Lets define helper function to plot 3D projection of the function def plot_trisurf(xfx, ax, legend="", elev=90, azim=-100, **kwargs): from matplotlib import cm """ Helper function to plot a 3D surface using a trisurf plot. Parameters: - xfx: A tuple containing the x-coordinates (2D points) and their corresponding function values. - ax: The matplotlib axis object for plotting. - legend: The legend/title for the plot. - elev, azim: Elevation and azimuth angles for the 3D view. - kwargs: Additional keyword arguments for further customization. """ xp, fxp = xfx[0], xfx[1] x, fx = xp, fxp X, Y = x[:, 0], x[:, 1] Z = fx.flatten() ax.plot_trisurf(X, Y, Z, antialiased=False, cmap=cm.jet) ax.view_init(azim=azim, elev=elev) ax.title.set_text(legend) # import CodPy's core module and Kernel class from codpy import core from codpy.kernel import Kernel .. GENERATED FROM PYTHON SOURCE LINES 173-183 Inverse Laplace operator In the experiments below, we will compare the original function with the result of applying the inverse Laplace operator to the function and then applying the Laplace operator again. This should yield the original function back, demonstrating the operator's effectiveness: $$ \\Delta_k^{-1}(\\Delta_k(f)) = f $$ ######################################################################## .. GENERATED FROM PYTHON SOURCE LINES 183-220 .. code-block:: Python def fun_Delta1(size_x=50, size_y=50): x, fx, z, fz, _, nabla_fz = generate_periodic_data_cartesian( size_x, size_y, periodic_fun, nabla_fun=nabla_my_fun ) nabla_fz = nabla_fz.reshape(-1, 2, 1) kernel_ptr = Kernel( x=x, fx=fx, set_kernel=core.kernel_setter("tensornorm", "scale_to_unitcube") ).get_kernel() temp = core.DiffOps.nabla_t_nabla_inv( x=x, y=x, fx=core.DiffOps.nabla_t_nabla(x=x, y=x, fx=fx, kernel_ptr=kernel_ptr), kernel_ptr=kernel_ptr, order=2, regularization=1e-8, ) multi_plot( [(x, fx), (x, temp)], plot_trisurf, projection="3d", mp_nrows=1, mp_figsize=(12, 3), mp_title=[ "Comparison between original function to the product of Laplace and its inverse" ], ) plt.show() fun_Delta1() .. image-sg:: /auto_ch3/images/sphx_glr_ch3_4_3_001.png :alt: ['Comparison between original function to the product of Laplace and its inverse'] :srcset: /auto_ch3/images/sphx_glr_ch3_4_3_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 221-233 Inverse Laplace operator In the experiments below, we will compare the original function with the result of applying the inverse Laplace operator to the function and then applying the Laplace operator again. This should yield the original function back, demonstrating the operator's effectiveness: $$ \\Delta_k(\\Delta_k^{-1}(f)) = f $$ Since the Laplace operator is associative, they should be equal. ######################################################################## .. GENERATED FROM PYTHON SOURCE LINES 233-266 .. code-block:: Python def fun_Delta2(size_x=50, size_y=50): x, fx, z, fz, _, nabla_fz = generate_periodic_data_cartesian( size_x, size_y, periodic_fun, nabla_fun=nabla_my_fun ) nabla_fz = nabla_fz.reshape(-1, 2, 1) kernel_ptr = Kernel( x=x, fx=fx, set_kernel=core.kernel_setter("tensornorm", "scale_to_unitcube") ).get_kernel() temp = core.DiffOps.nabla_t_nabla_inv( x=x, y=x, fx=fx, kernel_ptr=kernel_ptr, order=2, regularization=1e-8 ) temp = core.DiffOps.nabla_t_nabla( x=x, y=x, fx=temp, kernel_ptr=kernel_ptr, order=2, regularization=1e-8 ) multi_plot( [(x, fx), (x, temp)], plot_trisurf, projection="3d", mp_nrows=1, mp_figsize=(12, 3), mp_title=[ "Comparison between original function and the product of the inverse of the Laplace operator and the Laplace operator" ], ) plt.show() fun_Delta2() .. image-sg:: /auto_ch3/images/sphx_glr_ch3_4_3_002.png :alt: ['Comparison between original function and the product of the inverse of the Laplace operator and the Laplace operator'] :srcset: /auto_ch3/images/sphx_glr_ch3_4_3_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.527 seconds) .. _sphx_glr_download_auto_ch3_ch3_4_3.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ch3_4_3.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ch3_4_3.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ch3_4_3.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_