![](\"https://raw.githack.com/spectralDNS/spectralutilities/master/movies/Kuramato_movie_128.gif\")
Figure 1: Movie showing the evolution of the solution of the Kuramato-Sivashinsky equation.
\n", "" ] }, { "cell_type": "markdown", "id": "1b82a490", "metadata": { "editable": true }, "source": [ "## The Kuramato-Sivashinsky equation\n", "\n", "The Kuramato-Sivashinsky (KS) equation is known for its chaotic bahaviour, and it is\n", "often used in study of turbulence or turbulent combustion. We will here solve\n", "the KS equation in a doubly periodic domain $\\Omega=[-30\\pi, 30\\pi)^2$, starting from a\n", "single Gaussian pulse" ] }, { "cell_type": "markdown", "id": "40f1ebeb", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u(\\boldsymbol{x},t)}{\\partial t} + \\nabla^2 u(\\boldsymbol{x},t) + \\nabla^4\n", "u(\\boldsymbol{x},t) + |\\nabla u(\\boldsymbol{x},t)|^2 = 0 \\quad \\text{for }\\, \\boldsymbol{x} \\in \\Omega\n", "\\label{eq:ks} \\tag{1} \n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "0bf3843e", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u(\\boldsymbol{x}, 0) = \\exp(-0.01 \\boldsymbol{x} \\cdot \\boldsymbol{x}) \\notag\n", "\\label{_auto1} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "a3f91340", "metadata": { "editable": true }, "source": [ "## Spectral Galerkin method\n", "\n", "\n", "The PDE in ([1](#eq:ks)) can be solved with many different\n", "numerical methods. We will here use the [shenfun](https://github.com/spectralDNS/shenfun) software and this software makes use of\n", "the spectral Galerkin method. Being a Galerkin method, we need to reshape the\n", "governing equations into proper variational forms, and this is done by\n", "multiplying ([1](#eq:ks)) with the complex conjugate of a proper\n", "test function and then integrating\n", "over the domain. To this end we use testfunction $v\\in W^N(\\Omega)$, where $W^N(\\Omega)$\n", "is a suitable function space (defined in Eq. ([7](#eq:Wn))), and obtain" ] }, { "cell_type": "markdown", "id": "36175e5e", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial}{\\partial t} \\int_{\\Omega} u\\, \\overline{v}\\, w \\,dx = -\\int_{\\Omega}\n", "\\left(\\nabla^2 u + \\nabla^4 u \\ + |\\nabla u|^2 \\right) \\overline{v} \\, w \\,dx.\n", "\\label{eq:du_var} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "0b70e9df", "metadata": { "editable": true }, "source": [ "Note that the overline is used to indicate a complex conjugate, whereas $w$\n", "is a weight function. The function $u$\n", "is now to be considered a trial function, and the integrals over the\n", "domain are often referred to as inner products. With inner product notation" ] }, { "cell_type": "markdown", "id": "ce485733", "metadata": { "editable": true }, "source": [ "$$\n", "\\left(u, v\\right) = \\int_{\\Omega} u \\, \\overline{v} \\, w \\, dx.\n", "$$" ] }, { "cell_type": "markdown", "id": "571d5c61", "metadata": { "editable": true }, "source": [ "the variational problem can be formulated as" ] }, { "cell_type": "markdown", "id": "e0115f99", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial}{\\partial t} (u, v) = -\\left(\\nabla^2 u + \\nabla^4 u + |\\nabla u|^2,\n", "v \\right). \\label{eq:du_var2} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "4bd736f6", "metadata": { "editable": true }, "source": [ "The space and time discretizations are\n", "still left open. There are numerous different approaches that one could take for\n", "discretizing in time. Here we will use a fourth order exponential Runge-Kutta\n", "method." ] }, { "cell_type": "markdown", "id": "a3163d74", "metadata": { "editable": true }, "source": [ "## Discretization\n", "\n", "We discretize the model equation in space using continuously differentiable\n", "Fourier basis functions" ] }, { "cell_type": "markdown", "id": "1316d7f3", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\phi_l(x) = e^{\\imath \\underline{l} x}, \\quad -\\infty < l < \\infty,\n", "\\label{_auto2} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "cbfa19f1", "metadata": { "editable": true }, "source": [ "where $l$ is the wavenumber, and $\\underline{l}=\\frac{2\\pi}{L}l$ is the scaled wavenumber, scaled with domain\n", "length $L$ (here $60\\pi$). Since we want to solve these equations on a computer, we need to choose\n", "a finite number of test functions. A discrete function space $V^N$ can be defined as" ] }, { "cell_type": "markdown", "id": "fdce7ce9", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "V^N(x) = \\text{span} \\{\\phi_l(x)\\}_{l\\in \\boldsymbol{l}}, \\label{eq:Vn} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "621c6e53", "metadata": { "editable": true }, "source": [ "where $N$ is chosen as an even positive integer and $\\boldsymbol{l} = (-N/2,\n", "-N/2+1, \\ldots, N/2-1)$. And now, since $\\Omega$ is a\n", "two-dimensional domain, we can create a tensor product of two such one-dimensional\n", "spaces:" ] }, { "cell_type": "markdown", "id": "d595ac98", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "W^{\\boldsymbol{N}}(x, y) = V^N(x) \\otimes V^N(y), \\label{eq:Wn} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "aa8262e9", "metadata": { "editable": true }, "source": [ "where $\\boldsymbol{N} = (N, N)$. Obviously, it is not necessary to use the\n", "same number ($N$) of basis functions for each direction, but it is done here\n", "for simplicity. A 2D tensor product basis function is now defined as" ] }, { "cell_type": "markdown", "id": "46887195", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\Phi_{lm}(x,y) = e^{\\imath \\underline{l} x} e^{\\imath \\underline{m} y}\n", "= e^{\\imath (\\underline{l}x + \\underline{m}y )},\n", "\\label{_auto3} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "90815bf3", "metadata": { "editable": true }, "source": [ "where the indices for $y$-direction are $\\underline{m}=\\frac{2\\pi}{L}m$, and\n", "$\\boldsymbol{m}$ is the same set as $\\boldsymbol{l}$ due to using the same number of basis functions for each direction. One\n", "distinction, though, is that for the $y$-direction expansion coefficients are only stored for\n", "$m=(0, 1, \\ldots, N/2)$ due to Hermitian symmetry (real input data).\n", "\n", "We now look for solutions of the form" ] }, { "cell_type": "markdown", "id": "74359e27", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(x, y) = \\sum_{l=-N/2}^{N/2-1}\\sum_{m=-N/2}^{N/2-1}\n", "\\hat{u}_{lm} \\Phi_{lm}(x,y).\n", "\\label{_auto4} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "4d8a5f2b", "metadata": { "editable": true }, "source": [ "The expansion coefficients $\\hat{u}_{lm}$ can be related directly to the solution $u(x,\n", "y)$ using Fast Fourier Transforms (FFTs) if we are satisfied with obtaining\n", "the solution in quadrature points corresponding to" ] }, { "cell_type": "markdown", "id": "fa50c86c", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", " x_i = \\frac{60 \\pi i}{N}-30\\pi \\quad \\forall \\, i \\in \\boldsymbol{i},\n", "\\text{where}\\, \\boldsymbol{i}=(0,1,\\ldots,N-1), \n", "\\label{_auto5} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "66c43e0d", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", " y_j = \\frac{60 \\pi j}{N}-30\\pi \\quad \\forall \\, j \\in \\boldsymbol{j},\n", "\\text{where}\\, \\boldsymbol{j}=(0,1,\\ldots,N-1).\n", "\\label{_auto6} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "bdd44fd7", "metadata": { "editable": true }, "source": [ "Note that these points are different from the standard (like $2\\pi j/N$) since\n", "the domain\n", "is set to $[-30\\pi, 30\\pi]^2$ and not the more common $[0, 2\\pi]^2$. We now have" ] }, { "cell_type": "markdown", "id": "ae46333a", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\boldsymbol{u} =\n", "\\mathcal{F}_x^{-1}\\left(\\mathcal{F}_y^{-1}\\left(\\boldsymbol{\\hat{u}}\\right)\\right),\n", "\\label{_auto7} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "9ce62293", "metadata": { "editable": true }, "source": [ "where $\\boldsymbol{u} = \\{u(x_i, y_j)\\}_{(i,j)\\in \\boldsymbol{i} \\times \\boldsymbol{j}}$,\n", "$\\boldsymbol{\\hat{u}} = \\{\\hat{u}_{lm}\\}_{(l,m)\\in \\boldsymbol{l} \\times \\boldsymbol{m}}$\n", "and $\\mathcal{F}_x^{-1}$ is the inverse Fourier transform along direction\n", "$x$, for all indices in the other direction. Note that the two\n", "inverse FFTs are performed sequentially, one direction at the time, and that\n", "there is no scaling factor due\n", "the definition used for the inverse\n", "[Fourier transform](https://mpi4py-fft.readthedocs.io/en/latest/dft.html):" ] }, { "cell_type": "markdown", "id": "2da0c841", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u(x_j) = \\sum_{l=-N/2}^{N/2-1} \\hat{u}_l e^{\\imath \\underline{l}\n", "x_j}, \\quad \\,\\, \\forall \\, j \\in \\, \\boldsymbol{j}.\n", "\\label{_auto8} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "bde1ab80", "metadata": { "editable": true }, "source": [ "Note that this differs from the definition used by, e.g.,\n", "[Numpy](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html).\n", "\n", "The inner products used in Eq. ([4](#eq:du_var2)) may be\n", "computed using forward FFTs (using weight functions $w=1/L$):" ] }, { "cell_type": "markdown", "id": "ff7d3819", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\boldsymbol{\\hat{u}} =\n", "\\frac{1}{N^2}\n", "\\mathcal{F}_y\\left(\\mathcal{F}_x\\left(\\boldsymbol{u}\\right)\\right),\n", "\\label{_auto9} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "e3dcec3c", "metadata": { "editable": true }, "source": [ "From this we see that the variational forms\n", "may be written in terms of the Fourier transformed $\\hat{u}$. Expanding the\n", "exact derivatives of the nabla operator, we have" ] }, { "cell_type": "markdown", "id": "125917d9", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "(\\nabla^2 u, v) =\n", "\\left(-(\\underline{l}^2+\\underline{m}^2)\\hat{u}_{lm}\\right), \n", "\\label{_auto10} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "de4984ae", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "(\\nabla^4 u, v) = \\left((\\underline{l}^2+\\underline{m}^2)^2\\hat{u}_{lm}\\right), \n", "\\label{_auto11} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "c7b90c06", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "(|\\nabla u|^2, v) = \\left(\\widehat{|\\nabla u|^2}_{lm}\\right),\n", "\\label{_auto12} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "72605910", "metadata": { "editable": true }, "source": [ "where the indices on the right hand side run over $\\boldsymbol{l} \\times \\boldsymbol{m}$.\n", "We find that the equation to be solved for each wavenumber can be found directly as" ] }, { "cell_type": "markdown", "id": "f50b3d25", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial \\hat{u}_{lm}}{\\partial t} =\n", "\\left(\\underline{l}^2+\\underline{m}^2 -\n", "(\\underline{l}^2+\\underline{m}^2)^2\\right)\\hat{u}_{lm} - \\widehat{|\\nabla u|^2}_{lm},\n", "\\label{eq:du_var3} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "20229c39", "metadata": { "editable": true }, "source": [ "## Implementation\n", "\n", "The model equation ([1](#eq:ks)) is implemented in shenfun using Fourier basis functions for\n", "both $x$ and $y$ directions. We start the solver by implementing necessary\n", "functionality from required modules like [Numpy](https://numpy.org), [Sympy](https://sympy.org)\n", "and [matplotlib](https://matplotlib.org), in\n", "addition to [shenfun](https://github.com/spectralDNS/shenfun):" ] }, { "cell_type": "code", "execution_count": 1, "id": "0c1a7c93", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:01.718291Z", "iopub.status.busy": "2026-06-03T08:50:01.717983Z", "iopub.status.idle": "2026-06-03T08:50:02.256204Z", "shell.execute_reply": "2026-06-03T08:50:02.255815Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "from sympy import symbols, exp, lambdify\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import time\n", "from shenfun import *" ] }, { "cell_type": "markdown", "id": "b9d76d97", "metadata": { "editable": true }, "source": [ "The size of the problem (in real space) is then specified, before creating\n", "the [TensorProductSpace](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.tensorproductspace.TensorProductSpace), which is using a tensor product of two\n", "one-dimensional Fourier function spaces. We also\n", "create a [VectorSpace](https://shenfun.readthedocs.io/en/latest/shenfun.html#shenfun.tensorproductspace.VectorSpace), since this is required for computing the\n", "gradient of the scalar field `u`. The gradient is required for the nonlinear\n", "term." ] }, { "cell_type": "code", "execution_count": 2, "id": "063c35c2", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:02.258056Z", "iopub.status.busy": "2026-06-03T08:50:02.257752Z", "iopub.status.idle": "2026-06-03T08:50:02.780748Z", "shell.execute_reply": "2026-06-03T08:50:02.780486Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "# Size of discretization\n", "N = (128, 128)\n", "\n", "K0 = FunctionSpace(N[0], 'F', domain=(-30*np.pi, 30*np.pi), dtype='D')\n", "K1 = FunctionSpace(N[1], 'F', domain=(-30*np.pi, 30*np.pi), dtype='d')\n", "T = TensorProductSpace(comm, (K0, K1), **{'planner_effort': 'FFTW_MEASURE'})\n", "TV = VectorSpace([T, T])\n", "Tp = T.get_dealiased((1.5, 1.5))\n", "TVp = VectorSpace(Tp)" ] }, { "cell_type": "markdown", "id": "a3411ce3", "metadata": { "editable": true }, "source": [ "Test and trialfunctions are required for assembling the variational forms:" ] }, { "cell_type": "code", "execution_count": 3, "id": "3102955a", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:02.782276Z", "iopub.status.busy": "2026-06-03T08:50:02.782197Z", "iopub.status.idle": "2026-06-03T08:50:02.783877Z", "shell.execute_reply": "2026-06-03T08:50:02.783656Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "u = TrialFunction(T)\n", "v = TestFunction(T)" ] }, { "cell_type": "markdown", "id": "7ad972cb", "metadata": { "editable": true }, "source": [ "and some arrays are required to hold the solution. We also create an array\n", "`gradu`, that will be used to compute the gradient in the nonlinear term.\n", "Finally, the wavenumbers are collected in an array `K`. Here one feature is worth\n", "mentioning. The gradient in spectral space can be computed as `1j*K*U_hat`.\n", "However, since this is an odd derivative, and we are using an even number `N`\n", "for the size of the domain, the highest wavenumber must be set to zero. This is\n", "the purpose of the last keyword argument to `local_wavenumbers` below." ] }, { "cell_type": "code", "execution_count": 4, "id": "b9035407", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:02.785118Z", "iopub.status.busy": "2026-06-03T08:50:02.785050Z", "iopub.status.idle": "2026-06-03T08:50:02.796348Z", "shell.execute_reply": "2026-06-03T08:50:02.796101Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "x, y = symbols(\"x,y\", real=True)\n", "ue = exp(-0.01*(x**2+y**2))\n", "U = Array(T, buffer=ue)\n", "U_hat = Function(T)\n", "U_hat = U.forward(U_hat)\n", "mask = T.get_mask_nyquist()\n", "U_hat.mask_nyquist(mask)\n", "gradu = Array(TVp)\n", "K = np.array(T.local_wavenumbers(True, True, eliminate_highest_freq=True))\n", "X = T.local_mesh(True)" ] }, { "cell_type": "markdown", "id": "5681887e", "metadata": { "editable": true }, "source": [ "Note that using this `K` in computing derivatives has the same effect as\n", "achieved by symmetrizing the Fourier series by replacing the first sum below\n", "with the second when computing odd derivatives." ] }, { "cell_type": "markdown", "id": "573e1d22", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation}\n", "u = \\sum_{k=-N/2}^{N/2-1} \\hat{u}_k e^{\\imath k x}\n", "\\label{_auto13} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "434519b9", "metadata": { "editable": true }, "source": [ "\n", "\n", "\n", "$$\n", "\\begin{equation} \n", "u = \\sideset{}{'}\\sum_{k=-N/2}^{N/2} \\hat{u}_k e^{\\imath k x}\n", "\\label{_auto14} \\tag{20}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "id": "06e92d12", "metadata": { "editable": true }, "source": [ "Here $\\sideset{}{'}\\sum$ means that the first and last items in the sum are\n", "divided by two. Note that the two sums are equal as they stand (due to aliasing), but only the\n", "latter (known as the Fourier interpolant) gives the correct (zero) derivative of\n", "the basis with the highest wavenumber.\n", "\n", "Shenfun has a few integrators implemented in the [shenfun.utilities.integrators](https://shenfun.readthedocs.io/en/latest/shenfun.utilities.html#module-shenfun.utilities.integrators)\n", "submodule. Two such integrators are the 4th order explicit Runge-Kutta method\n", "`RK4`, and the exponential 4th order Runge-Kutta method `ETDRK4`. Both these\n", "integrators need two methods provided by the problem being solved, representing\n", "the linear and nonlinear terms in the problem equation. We define two methods\n", "below, called `LinearRHS` and `NonlinearRHS`" ] }, { "cell_type": "code", "execution_count": 5, "id": "80301cc6", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:02.798025Z", "iopub.status.busy": "2026-06-03T08:50:02.797936Z", "iopub.status.idle": "2026-06-03T08:50:02.800021Z", "shell.execute_reply": "2026-06-03T08:50:02.799809Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "def LinearRHS(self, u,**params):\n", " # Assemble diagonal bilinear forms\n", " return -(div(grad(u))+div(grad(div(grad(u)))))\n", "\n", "def NonlinearRHS(self, U, U_hat, dU, gradu, **params):\n", " # Assemble nonlinear term\n", " gradu = TVp.backward(1j*K*U_hat, gradu)\n", " dU = Tp.forward(0.5*(gradu[0]*gradu[0]+gradu[1]*gradu[1]), dU)\n", " dU.mask_nyquist(mask)\n", " dU *= -1\n", " return dU" ] }, { "cell_type": "markdown", "id": "762040be", "metadata": { "editable": true }, "source": [ "The code should, hopefully, be self-explanatory.\n", "\n", "All that remains now is to setup the\n", "integrator plus some plotting functionality for visualizing the results. Note\n", "that visualization is only nice when running the code in serial. For parallel,\n", "it is recommended to use [HDF5File](https://shenfun.readthedocs.io/en/latest/mpi4py_fft.io.html#mpi4py_fft.io.h5py_file.HDF5File), to store intermediate results to the HDF5\n", "format, for later viewing in, e.g., Paraview.\n", "\n", "We create an update function for plotting intermediate results with a\n", "cool colormap:" ] }, { "cell_type": "code", "execution_count": 6, "id": "8f63facd", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:02.801255Z", "iopub.status.busy": "2026-06-03T08:50:02.801179Z", "iopub.status.idle": "2026-06-03T08:50:02.803261Z", "shell.execute_reply": "2026-06-03T08:50:02.803046Z" }, "tags": [ "thebe-init" ] }, "outputs": [], "source": [ "from IPython.display import display\n", "\n", "# Integrate using an exponential time integrator\n", "def update(self, u, u_hat, t, tstep, **params):\n", " if tstep % params['plot_step'] == 0 and params['plot_step'] > 0:\n", " u = u_hat.backward(u)\n", " self.image = plt.contourf(X[0], X[1], u, 256, cmap=plt.get_cmap('hot'))\n", " self.image.axes.set_title(f'Energy {dx(u**2)}')\n", " display(self.image, clear=True)\n", " plt.pause(1e-6)" ] }, { "cell_type": "markdown", "id": "2d717922", "metadata": { "editable": true }, "source": [ "Now all that remains is to create the integrator and call it" ] }, { "cell_type": "code", "execution_count": 7, "id": "ab20ac9d", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2026-06-03T08:50:02.804505Z", "iopub.status.busy": "2026-06-03T08:50:02.804439Z", "iopub.status.idle": "2026-06-03T08:50:06.214862Z", "shell.execute_reply": "2026-06-03T08:50:06.214622Z" }, "tags": [ "thebe-init" ] }, "outputs": [ { "data": { "text/plain": [ "