{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Adaptive PDE discretizations on Cartesian grids\n", "## Volume : Divergence form PDEs\n", "## Part : One space dimension\n", "## Chapter : Heat and wave equations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We illustrate the discretization of time dependent partial differential equations in dimension one: diffusion (in divergence form) and the wave equation, whose PDE formulation read \n", "$$\n", " \\frac {\\partial u}{\\partial t} = \\frac {\\partial}{\\partial x} \\Big( c(x) \\frac {\\partial u}{\\partial x}\\Big),\n", " \\quad \\text{and} \\quad\n", " \\frac {\\partial^2 u}{\\partial t^2} = \\frac {\\partial}{\\partial x} \\Big( c(x) \\frac {\\partial u}{\\partial x}\\Big).\n", "$$\n", "For simplicity, we use periodic boundary conditions for the space variable, on the interval $\\Omega = [0,1]$. The diffusion coefficient $c : \\Omega \\to ]0,\\infty[$ is continuous and positive (it may also depend on time). \n", "\n", "We present the two most classical discretizations of the heat equation, explicit and implicit, their variational interpretation and their stability analysis. We rely on the Hamiltonian formalism for the wave equation.\n", "\n", "**Disclaimer.** This notebook does *not* contain original research. It is limited to elementary examples, and may serve as a gentle introduction to (some of) the numerical tools and techniques related to *time discretization*.\n", "My research is mainly devoted to the *spatial discretization* of PDE operators, which is irrelevant here in dimension one. Examples in dimension two and higher, involving non-trivial geometrical constructions, can be easily produced by combining the contents of the other notebooks with this one, and may also be presented in subsequent notebooks.\n", "\n", "**Stability analysis** We chose to present a stability analysis of the numerical schemes based on energetic arguments ($H^1$ semi-norm). There exists an alternative approach, perhaps slightly simpler, based on the Fourier transform. However, it requires the PDE coefficients to be constant, which is a strong limitation and defeats the purpose of this series of notebooks (discretizing PDEs with strong an position dependent anisotropy).\n", "\n", "**Related.** The discretization of time dependent *non-divergence* form PDEs is discussed [here](../Notebooks_NonDiv/Time1D_NonDiv.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[**Summary**](Summary.ipynb) of volume Divergence form PDEs, this series of notebooks.\n", "\n", "[**Main summary**](../Summary.ipynb) of the Adaptive Grid Discretizations \n", "\tbook of notebooks, including the other volumes.\n", "\n", "# Table of contents\n", " * [1. Quadratic forms and their discretization](#1.-Quadratic-forms-and-their-discretization)\n", " * [1.1 Continuous setting](#1.1-Continuous-setting)\n", " * [1.2 Discrete setting](#1.2-Discrete-setting)\n", " * [2. The heat equation, explicit scheme](#2.-The-heat-equation,-explicit-scheme)\n", " * [2.1 Discretization using automatic differentiation](#2.1-Discretization-using-automatic-differentiation)\n", " * [2.2 Stability analysis](#2.2-Stability-analysis)\n", " * [2.3 Sparse matrix](#2.3-Sparse-matrix)\n", " * [3. The heat equation, implicit scheme](#3.-The-heat-equation,-implicit-scheme)\n", " * [3.1 Discretization using automatic differentiation](#3.1-Discretization-using-automatic-differentiation)\n", " * [3.2 Stability analysis](#3.2-Stability-analysis)\n", " * [3.3 Sparse matrix](#3.3-Sparse-matrix)\n", " * [4. Wave equation](#4.-Wave-equation)\n", " * [4.1 Discretization using automatic differentiation](#4.1-Discretization-using-automatic-differentiation)\n", " * [4.2 Stability analysis](#4.2-Stability-analysis)\n", " * [4.3 Sparse matrices](#4.3-Sparse-matrices)\n", " * [4.4 Using the Hamiltonian class](#4.4-Using-the-Hamiltonian-class)\n", "\n", "\n", "\n", "**Acknowledgement.** Some of the experiments presented in these notebooks are part of \n", "ongoing research with Ludovic Métivier and Da Chen.\n", "\n", "Copyright Jean-Marie Mirebeau, Centre Borelli, ENS Paris-Saclay, CNRS, University Paris-Saclay" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 0. Importing the required libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.434529Z", "iopub.status.busy": "2024-04-30T09:15:13.434263Z", "iopub.status.idle": "2024-04-30T09:15:13.439478Z", "shell.execute_reply": "2024-04-30T09:15:13.439067Z" } }, "outputs": [], "source": [ "import sys; sys.path.insert(0,\"..\") # Allow import of agd from parent directory (useless if conda package installed)\n", "#from Miscellaneous import TocTools; print(TocTools.displayTOC('Time1D_Div','Div'))" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.441742Z", "iopub.status.busy": "2024-04-30T09:15:13.441603Z", "iopub.status.idle": "2024-04-30T09:15:13.720743Z", "shell.execute_reply": "2024-04-30T09:15:13.720282Z" } }, "outputs": [], "source": [ "from agd import FiniteDifferences as fd\n", "from agd import AutomaticDifferentiation as ad\n", "from agd.Plotting import animation_curve\n", "from agd.ODE.hamiltonian import QuadraticHamiltonian\n", "norm_infinity = ad.Optimization.norm_infinity " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.723311Z", "iopub.status.busy": "2024-04-30T09:15:13.723051Z", "iopub.status.idle": "2024-04-30T09:15:13.725847Z", "shell.execute_reply": "2024-04-30T09:15:13.725348Z" } }, "outputs": [], "source": [ "import numpy as np\n", "import scipy.sparse\n", "import matplotlib.pyplot as plt\n", "from matplotlib import rc; rc('animation', html='html5')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some utility functions" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.727979Z", "iopub.status.busy": "2024-04-30T09:15:13.727734Z", "iopub.status.idle": "2024-04-30T09:15:13.730786Z", "shell.execute_reply": "2024-04-30T09:15:13.730406Z" } }, "outputs": [], "source": [ "from agd.ExportedCode.Notebooks_NonDiv.Time1D_NonDiv import accumulate # Or itertools.accumulate if python>=3.8" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.732384Z", "iopub.status.busy": "2024-04-30T09:15:13.732265Z", "iopub.status.idle": "2024-04-30T09:15:13.734354Z", "shell.execute_reply": "2024-04-30T09:15:13.734089Z" } }, "outputs": [], "source": [ "def reload_packages():\n", " from Miscellaneous.rreload import rreload\n", " global fd,ad,QuadraticHamiltonian\n", " fd,ad,QuadraticHamiltonian = rreload([fd,ad,QuadraticHamiltonian],rootdir=\"..\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Quadratic forms and their discretization\n", "\n", "We introduce the basic objects underlying the heat and wave equations, which are quadratic forms on suitable spaces of functions. We then introduce their discretization." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Continuous setting\n", "\n", "We briefly recall some of the concepts underlying divergence form PDE discretizations. The objective here is to set notations, not present a course on the subject, and not to achieve absolute mathematical rigor either. Please consider a textbook for that purpose.\n", "\n", "**Quadratic forms, polarization and associated operator.**\n", "A bilinear symmetric form is a function $Q(u,v)$ of all elements $u,v$ of a vector space $V$, which is linear w.r.t \n", "$u$ and $v$ (separately, not simultaneously), and obeys the symmetry $Q(u,v)=Q(v,u)$ for all $u,v \\in V$.\n", "\n", "A quadratic form is the specialization $Q(u,u)$ of a bilinear symmetric form to the case $u=v$, $u\\in V$. No information is lost in this process, since the bilinear symmetric form can be recovered using an identity referred to as polarization. The quadratic form is said non-negative if $Q(u,u)\\geq 0$ for all $u\\in V$.\n", "\n", "One can associate a linear operator to a quadratic form $Q$ whose domain $V$ is a (dense subset of a) Hilbert space $H$. This operator is denoted by the same letter, and defined by the identity\n", "$$\n", " = Q(u,v)\n", "$$\n", "for all $u,v \\in V$. \n", "\n", "**Case of the $H^1$ semi-norm.** We denote by $Q$ the Dirichlet elliptic energy, defined for all $u \\in H^1(\\Omega)$ by \n", "$$\n", " Q(u,u) := \\int_\\Omega c(x)|\\nabla u(x)|^2 \\ dx.\n", "$$\n", "Note that $Q$ is non-negative, since $c\\geq 0$.\n", "The ambient space is $L^2(\\Omega)$, equipped with the scalar product defined by \n", "$$\n", " = \\int_\\Omega u(x)v(x) \\ dx.\n", "$$\n", "\n", "**PDE reformulation.**\n", "An integration by parts shows that, recalling that $\\Omega=[0,1]$ with periodic b.c.,\n", "$$\n", " Q u = - \\frac {\\partial}{\\partial x} \\Big( c(x) \\frac {\\partial u}{\\partial x}\\Big)\n", "$$\n", "Therefore the heat and wave equations can be reformulated as \n", "$$\n", " \\frac d {dt} u = - Q u,\n", " \\quad \\text{and} \\quad\n", " \\frac {d^2} {dt^2} u = - Q u. \n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Discrete setting\n", "\n", "For any $h>0$, which is the inverse of a positive integer, we introduce the set $\\Omega_h := \\Omega\\cap hZ$. In the numerical codes, $h$ is denoted `dx`.\n", "\n", "**Discretization.**\n", "Define the quadratic form $Q_h$ by \n", "$$\n", " Q_h(u,u) := \\sum_{x \\in \\Omega_h} \\frac{c(x)} 2\\Big[\\big(\\frac {u(x+h)-u(x)} h\\big)^2+\\big(\\frac {u(x-h)-u(x)} h\\big)^2\\Big]\n", "$$\n", "as well as $I_h$ defined by \n", "$$\n", " I_h(u,u) := \\sum_{x \\in \\Omega_h} u(x)^2.\n", "$$\n", "\n", "**Consistency.** For twice continuously differentiable $u,v:\\Omega \\to R$, one has second order consistency\n", "$$\n", " h\\ Q_h(u,v) = Q(u,v) + O(h^2),\n", " \\quad \\text{and} \\quad\n", " h\\ I_h(u,v) = + O(h^2).\n", "$$\n", "\n", "We could have included the factor $h$ in the definition of $Q_h$ and $I_h$, but chose not to for simplicity. (By eliminating this factor, the matrix of $I_h$ is the identity, and we recover the usual scalar product on $R^d$. \n", "Note that this simplification would not be possible in the finite element setting, where $I_h$ is known as the mass matrix.\n", "In addition, this factor would need to be adjusted to $h^d$ in dimension $d\\geq 1$.)\n", "\n", "Note finally that the gradient, and gradient flow, are not modified, if the energy functional and the scalar product are multiplied by an *identical* positive constant, here the constant $h>0$.\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us introduce the discretized operators $Q_h$ and $I_h$.\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.736175Z", "iopub.status.busy": "2024-04-30T09:15:13.736066Z", "iopub.status.idle": "2024-04-30T09:15:13.739310Z", "shell.execute_reply": "2024-04-30T09:15:13.738974Z" } }, "outputs": [], "source": [ "def I(u,v):\n", " \"\"\"Approximation of the L2 scalar product\"\"\"\n", " return np.sum(u*v,axis=0)\n", "\n", "def Q(u,v,c,dx):\n", " \"\"\"\n", " Finite differences discretization of the H1 bilinear form,\n", " in dimension 1, with periodic b.c.\n", " \"\"\"\n", " dup = fd.DiffUpwind(u,( 1,),dx,padding=None) # (u(x+dx)-u(x))/dx\n", " dvp = fd.DiffUpwind(v,( 1,),dx,padding=None) # (v(x+dx)-v(x))/dx\n", "\n", " dum = fd.DiffUpwind(u,(-1,),dx,padding=None) # (u(x-dx)-u(x))/dx\n", " dvm = fd.DiffUpwind(v,(-1,),dx,padding=None) # (v(x-dx)-u(x))/dx\n", " \n", " r = 0.5*c*(dup*dvp+dum*dvm)\n", " return r.sum(axis=0) # Integate r on domain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For concreteness, we introduce a discretization grid, some diffusion coefficients, and test functions." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.741328Z", "iopub.status.busy": "2024-04-30T09:15:13.741185Z", "iopub.status.idle": "2024-04-30T09:15:13.743444Z", "shell.execute_reply": "2024-04-30T09:15:13.743145Z" } }, "outputs": [], "source": [ "X,dx = np.linspace(0,1,100,endpoint=False,retstep=True)\n", "Tmax = 0.5 # Default time interval is [0,Tmax]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.744957Z", "iopub.status.busy": "2024-04-30T09:15:13.744843Z", "iopub.status.idle": "2024-04-30T09:15:13.747618Z", "shell.execute_reply": "2024-04-30T09:15:13.747161Z" } }, "outputs": [], "source": [ "pi2 = np.pi*2.\n", "c_constant = 0.7\n", "c_positive = 0.7 + 0.4*np.sin(pi2*X)\n", "\n", "u_disc = 1.*(X>=0.5)*(X<=0.75)\n", "u_cont = np.maximum(0.,(0.5-X)*(X-0.75)); u_cont/=np.max(u_cont)\n", "u_smooth = u_cont**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can test our bilinear forms on arbirary numpy arrays." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.749970Z", "iopub.status.busy": "2024-04-30T09:15:13.749818Z", "iopub.status.idle": "2024-04-30T09:15:13.754188Z", "shell.execute_reply": "2024-04-30T09:15:13.753895Z" } }, "outputs": [ { "data": { "text/plain": [ "(-11.623155409290282, -0.42073549240394836)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u_ = np.sin(pi2*X)\n", "v_ = np.cos(pi2*X+1)\n", "dx*Q(u_,v_,c_positive,dx), dx*I(u_,v_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Specializing to $u=v$ yields a non-negative value, by construction." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.756051Z", "iopub.status.busy": "2024-04-30T09:15:13.755905Z", "iopub.status.idle": "2024-04-30T09:15:13.758983Z", "shell.execute_reply": "2024-04-30T09:15:13.758684Z" } }, "outputs": [ { "data": { "text/plain": [ "(13.812901002099066, 0.5)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dx*Q(u_,u_,c_positive,dx), dx*I(u_,u_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. The heat equation, explicit scheme\n", "\n", "\n", "**Continuous setting.** The heat equation can be written in the so-called variational form\n", "$$\n", " <\\frac d {dt} u, v> = - Q(u,v)\n", "$$\n", "for all test functions $v$. Recall that $<,>$ denotes the $L^2$ scalar product, and $Q$ the bilinear form associated with the Dirichlet energy.\n", "\n", "**Discretization.**\n", "We use a first order explicit time discretization, which yields the system\n", "$$\n", " I_h\\big(\\frac{u_{n+1}-u_n}{\\delta t},v\\big) = - Q_h(u_n,v)\n", "$$\n", "for all test functions $v$, where $u,v,w : \\Omega_h \\to R$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Discretization using automatic differentiation\n", "\n", "As a first approach, we rely on automatic differentiation to implement the variational formulation directly.\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.760850Z", "iopub.status.busy": "2024-04-30T09:15:13.760725Z", "iopub.status.idle": "2024-04-30T09:15:13.763468Z", "shell.execute_reply": "2024-04-30T09:15:13.763149Z" } }, "outputs": [], "source": [ "def HeatExplicit(u,c,dx,dt):\n", " \"\"\"One time step of the explicit scheme for the heat equation\"\"\"\n", " \n", " # Define independent second order AD variables\n", " u_ad = ad.Sparse2.identity(u.shape) # Unknown u_{n+1}\n", " v_ad = ad.Sparse2.identity(u.shape,shift=u.size) # Test function\n", " \n", " dtu_ad = (u_ad-u)/dt # Time derivative (u_{n+1}-u_n)/dt\n", " \n", " weakform = I(dtu_ad,v_ad) + Q(u,v_ad,c,dx) # Should vanish for all test functions v_ad\n", " \n", " return weakform.solve_weakform() # Finds u_{n+1}" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.765115Z", "iopub.status.busy": "2024-04-30T09:15:13.765010Z", "iopub.status.idle": "2024-04-30T09:15:13.834275Z", "shell.execute_reply": "2024-04-30T09:15:13.833814Z" } }, "outputs": [], "source": [ "dt = 4e-5\n", "solution =np.array(list(accumulate(\n", " np.arange(0,100*dt,dt), # time interval\n", " initial=u_disc, # initial condition\n", " func=lambda u,t: HeatExplicit(u,c_positive,dx,dt) # evolution rule\n", ")))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:13.836790Z", "iopub.status.busy": "2024-04-30T09:15:13.836658Z", "iopub.status.idle": "2024-04-30T09:15:14.999147Z", "shell.execute_reply": "2024-04-30T09:15:14.998801Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "animation_curve(X,solution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Stability analysis\n", "\n", "We can rephrase the update of the explicit scheme for the heat equation in the form\n", "$$\n", " u_{n+1} = (\\mathrm{Id} - I_h^{-1} Q_h \\delta t) u_n.\n", "$$\n", "where $\\mathrm{Id}$ denotes the identity matrix.\n", "\n", "The scheme is stable, in the $L^2$ norm, iff all the eigenvalues of $\\mathrm{Id} - I_h^{-1} Q_h \\delta t$ lie in $[-1,1]$. Equivalently iff, in the sense of symmetric matrices\n", "$$\n", "Q_h \\delta t \\preceq 2 I_h.\n", "$$\n", "Said otherwise, one must have $Q_h(u,u) \\delta t \\preceq 2 I_h(u,u)$ for any $u : \\Omega_h \\to R$. Based on the inequality $(a-b)^2 \\leq 2(a^2+b^2)$, and some simplifications, we obtain\n", "$$\n", " Q_h(u,u) \n", " \\leq h^{-1} \\sum_{x \\in \\Omega_h} c(x) \\big[u(x+h)^2+2 u(x)^2 +u(x-h)^2\\big] \n", " \\leq \\frac {4c_{\\max}} h \\sum_{x \\in \\Omega_h} u(x)^2 = \\frac {4c_{\\max}} {h^2} I_h(u,u),\n", "$$\n", "where $c_{\\max}$ is a uniform bound for the diffusion coefficient $c$. This results in the CFL condition\n", "$$\n", " 2c_{\\max}\\delta t \\leq h^2.\n", "$$\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.000785Z", "iopub.status.busy": "2024-04-30T09:15:15.000669Z", "iopub.status.idle": "2024-04-30T09:15:15.002579Z", "shell.execute_reply": "2024-04-30T09:15:15.002349Z" } }, "outputs": [], "source": [ "def HeatExplicit_CFL(c,dx):\n", " \"\"\"\n", " The explicit discretization of the heat equation\n", " is stable provided the time step dt is below this threshold.\n", " \"\"\"\n", " return dx**2/(2*np.max(c))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.003970Z", "iopub.status.busy": "2024-04-30T09:15:15.003873Z", "iopub.status.idle": "2024-04-30T09:15:15.005975Z", "shell.execute_reply": "2024-04-30T09:15:15.005726Z" } }, "outputs": [ { "data": { "text/plain": [ "4.545454545454545e-05" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "HeatExplicit_CFL(c_positive,dx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Sparse matrix \n", "\n", "The implementation presented above of the heat equation lacks in numerical efficiency, and perhaps explicitness, due to the heavy use of automatic differentiation.\n", "To remedy this problem, we extract the matrices involved, and replace the weak form solution with an explicit matrix-vector product." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us extract the (sparse symmetric) matrix associated with quadratic form $Q_h$.\n", "Note the $1/2$ factor, as the matrix of interest is *half the hessian* of $Q_h(u,u)$." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.007331Z", "iopub.status.busy": "2024-04-30T09:15:15.007235Z", "iopub.status.idle": "2024-04-30T09:15:15.010431Z", "shell.execute_reply": "2024-04-30T09:15:15.010140Z" } }, "outputs": [], "source": [ "u_ad = ad.Sparse2.identity(X.shape)\n", "Quu_ad = 0.5*Q(u_ad,u_ad,c_positive,dx) #Note 1/2 factor\n", "Q_mat = scipy.sparse.coo_matrix(Quu_ad.triplets()).tocsc() # Hessian matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The numerical solution is computed much faster. It is also identical up to floating point rounding error.\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.011758Z", "iopub.status.busy": "2024-04-30T09:15:15.011677Z", "iopub.status.idle": "2024-04-30T09:15:15.016009Z", "shell.execute_reply": "2024-04-30T09:15:15.015779Z" } }, "outputs": [], "source": [ "dt = 4e-5\n", "solution2 = np.array(list(accumulate(\n", " np.arange(0,100*dt,dt), # time interval\n", " initial=u_disc, # initial condition\n", " func=lambda u,t: u-dt*Q_mat*u # evolution rule\n", ")))" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.017327Z", "iopub.status.busy": "2024-04-30T09:15:15.017246Z", "iopub.status.idle": "2024-04-30T09:15:15.019072Z", "shell.execute_reply": "2024-04-30T09:15:15.018833Z" } }, "outputs": [], "source": [ "assert norm_infinity(solution-solution2) < 1e-14" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. The heat equation, implicit scheme\n", "\n", "**Continuous setting.**\n", "The heat equation can be regarded as the gradient flow of the energy $\\frac 1 2 Q(u,u)$ w.r.t. the $L^2$ metric. \n", "A semi-discretized formulation reflecting this principle is as follows: define $u_{n+1}$ as the minimizer to \n", "$$\n", " \\min_{u_{n+1}} \\frac 1 {\\delta t} \\|u_{n+1}-u_n\\|^2 + Q(u_{n+1},u_{n+1}).\n", "$$\n", "There exists another, non-linear, similarly looking formulation of the heat equation: as the gradient flow of the entropy w.r.t the Wasserstein metric.\n", "\n", "Each time step of the implicit scheme involves solving a linear equation, which is more expensive than the matrix-vector products involved in the explicit scheme. However, the time steps used in the implicit scheme are unconstrained, and the scheme is unconditionally stable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Discretization using automatic differentiation\n", "\n", "We translate in the discrete setting the variational principle defining $u_{n+1}$, in the form\n", "$$\n", " \\min_{u_{n+1}} \\frac 1 {\\delta t} I_h(u_{n+1}-u_n,\\ u_{n+1}-u_n) + \n", " Q_h(u_{n+1},u_{n+1}).\n", "$$\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.020575Z", "iopub.status.busy": "2024-04-30T09:15:15.020497Z", "iopub.status.idle": "2024-04-30T09:15:15.022458Z", "shell.execute_reply": "2024-04-30T09:15:15.022250Z" } }, "outputs": [], "source": [ "def HeatImplicit(u,c,dx,dt):\n", " u_ad = ad.Sparse2.identity(u.shape) # Unknown u_{n+1}\n", " energy = (1/dt)*I(u_ad-u,u_ad-u) + Q(u_ad,u_ad,c,dx)\n", " return energy.solve_stationnary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As announced, one of the key interests of the implicit scheme is the ability to take large time steps, here several orders of magnitude larger than for the explicit scheme." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.023865Z", "iopub.status.busy": "2024-04-30T09:15:15.023791Z", "iopub.status.idle": "2024-04-30T09:15:15.042773Z", "shell.execute_reply": "2024-04-30T09:15:15.042467Z" } }, "outputs": [], "source": [ "dt = 1e-2\n", "solution =np.array(list(accumulate(\n", " np.arange(0,20*dt,dt), # time interval\n", " initial=u_disc, # initial condition\n", " func=lambda u,t: HeatImplicit(u,c_positive,dx,dt) # evolution rule\n", ")))" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.044321Z", "iopub.status.busy": "2024-04-30T09:15:15.044231Z", "iopub.status.idle": "2024-04-30T09:15:15.301310Z", "shell.execute_reply": "2024-04-30T09:15:15.300976Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "animation_curve(X,solution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Stability analysis\n", "\n", "We can rephrase the update of the explicit scheme for the heat equation in the form\n", "$$\n", " u_{n+1} = (\\mathrm{Id} + I_h^{-1} Q_h \\delta t)^{-1} u_n.\n", "$$\n", "where $\\mathrm{Id}$ denotes the identity matrix. \n", "\n", "The scheme is stable provided all the eigenvalues of the operator $A = (\\mathrm{Id} + I_h^{-1} Q_h \\delta t)^{-1}$ lie in $[-1,1]$. This property holds unconditionally w.r.t $\\delta t \\geq 0$. Indeed, the eigenvalues of $A$ take the form\n", "$$\n", " \\frac 1 {1+\\lambda h/\\delta t}\n", "$$\n", "where $\\lambda$ is an eigenvalue of $Q_h$, hence $\\lambda \\geq 0$ since $Q_h$ is a non-negative symmetric matrix. (We use the fact that $I_h = h \\mathrm{Id}$.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Sparse matrix\n", "\n", "The implicit scheme for the heat equation involves solving repeatedly the same (sparse, positive definite) linear system. This procedure can be made numerically more efficient by pre-factorization of the matrix." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.302999Z", "iopub.status.busy": "2024-04-30T09:15:15.302875Z", "iopub.status.idle": "2024-04-30T09:15:15.305378Z", "shell.execute_reply": "2024-04-30T09:15:15.305118Z" } }, "outputs": [], "source": [ "impl_mat = scipy.sparse.eye(len(u_disc)) + dt*Q_mat\n", "solver = scipy.sparse.linalg.factorized(impl_mat.tocsc())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The scheme runs faster, but yields an identical solution up to machine precision." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.306795Z", "iopub.status.busy": "2024-04-30T09:15:15.306692Z", "iopub.status.idle": "2024-04-30T09:15:15.308780Z", "shell.execute_reply": "2024-04-30T09:15:15.308546Z" } }, "outputs": [], "source": [ "dt = 1e-2\n", "solution2 = np.array(list(accumulate(\n", " np.arange(0,20*dt,dt), # time interval\n", " initial=u_disc, # initial condition\n", " func=lambda u,t: solver(u) # evolution rule\n", ")))" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.310150Z", "iopub.status.busy": "2024-04-30T09:15:15.310056Z", "iopub.status.idle": "2024-04-30T09:15:15.311719Z", "shell.execute_reply": "2024-04-30T09:15:15.311493Z" } }, "outputs": [], "source": [ "assert norm_infinity(solution-solution2) < 1e-13" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Wave equation\n", "\n", "**Continuous setting.**\n", "The wave equation can be rephrased in the Hamiltonian formalism\n", "$$\n", " \\frac {d u} {dt} = \\frac{\\partial H}{\\partial \\dot u} \n", " \\quad \\text{and} \\quad \n", " \\frac {d\\dot u} {dt} = - \\frac{\\partial H}{\\partial u}.\n", "$$\n", "The Hamiltonian is defined as the sum of a *potential* elastic energy, and of a *kinetic* energy\n", "$$\n", " H(u,\\dot u) := \\frac 1 2 \\big( Q(u,u) + <\\dot u,\\dot u> \\big),\n", "$$\n", "where the bilinear form $Q$ and scalar product $<,>$ are those already introduced.\n", "The Hamitlonian is constant in time, along solutions of sufficient regularity.\n", "The discretization mimicks this formulation, using a symplectic integrator to conserve the Hamiltonian and ensure stability." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.1 Discretization using automatic differentiation\n", "\n", "**Translation of the continuous setting.**\n", "We introduce a discretized Hamiltonian, defined by\n", "$$\n", " H^0_h(u,\\dot u) := \\frac 1 2 \\big(Q_h(u,u) + I_h(\\dot u,\\dot u) \\big).\n", "$$\n", "The Euler symplectic integrator is defined as \n", "$$\n", " \\frac{u_{n+1}-u_n}{\\delta t} = \\frac{\\partial H^0_h}{\\partial \\dot u}(u_{n+1},\\dot u_n)\n", " \\quad \\text{and} \\quad\n", " \\frac{\\dot u_{n+1}-\\dot u_n}{\\delta t} = - \\frac{\\partial H^0_h}{\\partial u}(u_{n+1},\\dot u_n),\n", "$$\n", "where $u_n,\\dot u_n,u_{n+1},\\dot u_{n+1} : \\Omega_h \\to R$.\n", "In this special case, the variable $\\dot u_n$ does approximate the time derivative of $u$, hence the choice of notation. However, in general, when e.g. the medium density is not constant, the second variable in the Hamiltonian corresponds to the momentum.\n", "\n", "**Note on separability.**\n", "For a general Hamiltonian, the first step would be implicit, see for instance this [notebook](../Notebooks_Algo/Dense.ipynb). However, the wave equation Hamiltonian benefits from a separable structure, which makes the two updates explicit.\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.313094Z", "iopub.status.busy": "2024-04-30T09:15:15.313012Z", "iopub.status.idle": "2024-04-30T09:15:15.314741Z", "shell.execute_reply": "2024-04-30T09:15:15.314525Z" } }, "outputs": [], "source": [ "def WaveHamiltonian(u,up,c,dx):\n", " return 0.5 * ( Q(u,u,c,dx) + I(up,up) )" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.316075Z", "iopub.status.busy": "2024-04-30T09:15:15.315998Z", "iopub.status.idle": "2024-04-30T09:15:15.318140Z", "shell.execute_reply": "2024-04-30T09:15:15.317917Z" } }, "outputs": [], "source": [ "def SeparableSymplectic(q,p,H,dt,dx):\n", " \"\"\"Euler's symplectic integrator, for a Hamiltonian assumed to be separable\"\"\"\n", " # Using the separable structure here (q instead of q_next in r.h.s)\n", " p_ad = ad.Sparse.identity(constant=p)\n", " q_next = q+dt*H(q,p_ad).to_dense().gradient()\n", " \n", " q_ad = ad.Sparse.identity(constant=q_next)\n", " p_next = p-dt*H(q_ad,p).to_dense().gradient() \n", " return q_next,p_next" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the wave equation does not have a regularizing effect, unlike the heat equation, we use a continuous initial condition." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.319463Z", "iopub.status.busy": "2024-04-30T09:15:15.319385Z", "iopub.status.idle": "2024-04-30T09:15:15.433368Z", "shell.execute_reply": "2024-04-30T09:15:15.432986Z" } }, "outputs": [], "source": [ "def H(q,p): return WaveHamiltonian(q,p,c_positive,dx)\n", "dt = 9e-3\n", "solution = np.array(list(accumulate(\n", " np.arange(0,200*dt,dt), # time interval\n", " initial=(u_cont,0.*u_cont), # initial condition\n", " func=lambda state,t: SeparableSymplectic(*state,H,dt,dx) # state contains q and p\n", ")))" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:15.434933Z", "iopub.status.busy": "2024-04-30T09:15:15.434839Z", "iopub.status.idle": "2024-04-30T09:15:17.575203Z", "shell.execute_reply": "2024-04-30T09:15:17.574873Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "animation_curve(X,solution[:,0,:])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.2 Stability analysis\n", "\n", "**Separable quadratic Hamiltonians.**\n", "Symplectic integrators preserve a perturbation of the Hamiltonian, and this property is the key to stability analysis. We briefly present this theory in the case of separable quadratic Hamilonians on $R^N \\times R^N$, which is enough for our purposes. Assume that \n", "$$\n", " H(q,p) = \\frac 1 2 ( + ),\n", "$$\n", "where $A$ and $B$ are symmetric matrices. Euler's symplectic scheme yields\n", "$$\n", " q_{n+1} = q_n + \\delta t A p_n,\n", " \\quad \\text{and} \\quad\n", " p_{n+1} = p_n - \\delta t B q_n.\n", "$$\n", "An elementary computation shows that the following quantity is *exactly* conserved along the iterations\n", "$$\n", " \\tilde H(q,p) := \\frac 1 2 ( + + \\delta t )\n", "$$\n", "If $\\tilde H$ is a positive definite quadratic form, then the iterates $(q_n,p_n)_{n\\geq 0}$ are bounded independently of $n$, meaning that the scheme is stable. This condition amounts to \n", "$$\n", " \\begin{pmatrix}\n", " B & \\frac {\\delta t} 2 B A\\\\\n", " \\frac {\\delta t} 2 A B & A\n", " \\end{pmatrix}\n", " \\succ 0\n", "$$\n", "where $M\\succ N$ means that $M-N$ is positive definite. If $A$ and $B$ commute, then this condition can be simplified to \n", "$$\n", " B \\succ 0, \\quad \n", " A \\succ 0, \\quad \n", " \\delta t^2 BA \\prec 4 \\mathrm{Id}.\n", "$$\n", "\n", "**Specialization to the wave equation**\n", "We neglect here the technicalities associated with positive definiteness versus semi-definiteness.\n", "The operators $Q_h$ and $I_h$ are positive (semi-)definite and commute, thus only the third condition remains. After specialization we obtain the stability criterion\n", "$$\n", " \\delta t^2 Q_h \\preceq 4 I_h.\n", "$$\n", "(We used the identity $I_h = \\mathrm{Id}$.) From the same analysis as in the explicit Euler equation, based on the specific form of $Q_h$, we obtain the condition\n", "$$\n", " \\delta t^2 c_{\\max} \\leq \\delta x^2.\n", "$$\n", "(Note that the wave velocity is $\\sqrt{c_{\\max}}$.)\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.577523Z", "iopub.status.busy": "2024-04-30T09:15:17.577428Z", "iopub.status.idle": "2024-04-30T09:15:17.579419Z", "shell.execute_reply": "2024-04-30T09:15:17.579195Z" } }, "outputs": [], "source": [ "def Wave_CFL(c,dx):\n", " \"\"\"Returns the largest time step for which the wave equation scheme is stable.\"\"\"\n", " return dx/np.sqrt(np.max(c))" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.580751Z", "iopub.status.busy": "2024-04-30T09:15:17.580667Z", "iopub.status.idle": "2024-04-30T09:15:17.582882Z", "shell.execute_reply": "2024-04-30T09:15:17.582629Z" } }, "outputs": [ { "data": { "text/plain": [ "(0.009534625892455923, 0.009)" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Wave_CFL(c_positive,dx),dt" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.584177Z", "iopub.status.busy": "2024-04-30T09:15:17.584096Z", "iopub.status.idle": "2024-04-30T09:15:17.585947Z", "shell.execute_reply": "2024-04-30T09:15:17.585710Z" } }, "outputs": [], "source": [ "def ConservedHamiltonian(u,up,c,dx,dt):\n", " \"\"\"\n", " Perturbation of the discretized Wave Hamiltonian, \n", " which is exactly conserved by Euler's symplectic scheme.\n", " \"\"\"\n", " return WaveHamiltonian(u,up,c,dx) + dt*Q(u,up,c,dx)/2" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.587295Z", "iopub.status.busy": "2024-04-30T09:15:17.587219Z", "iopub.status.idle": "2024-04-30T09:15:17.642307Z", "shell.execute_reply": "2024-04-30T09:15:17.642064Z" } }, "outputs": [], "source": [ "invariant = [ConservedHamiltonian(u,up,c_positive,dx,dt) for u,up in solution]\n", "assert norm_infinity(np.array(invariant)-invariant[0])<1e-12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.3 Sparse matrices\n", "\n", "We improve the computational efficiency, and explicitness, of the numerical scheme by using sparse matrix-vector products instead of relying on automatic differentiation.\n", "\n", "\n", "\n", "\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.643753Z", "iopub.status.busy": "2024-04-30T09:15:17.643670Z", "iopub.status.idle": "2024-04-30T09:15:17.645511Z", "shell.execute_reply": "2024-04-30T09:15:17.645265Z" } }, "outputs": [], "source": [ "def WaveSymplectic(u,up,dt):\n", " \"\"\"One time step of the Euler symplectic scheme, updating u first\"\"\"\n", " u_next = u+dt*up\n", " up_next = up-dt*Q_mat*u_next # Q_mat is the matrix of Q_h with c_positive\n", " return u_next,up_next " ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.646904Z", "iopub.status.busy": "2024-04-30T09:15:17.646829Z", "iopub.status.idle": "2024-04-30T09:15:17.653836Z", "shell.execute_reply": "2024-04-30T09:15:17.653616Z" } }, "outputs": [], "source": [ "dt = 9e-3\n", "solution2 = np.array(list(accumulate(\n", " np.arange(0,200*dt,dt), # time interval\n", " initial=(u_cont,0.*u_cont), # initial condition\n", " func=lambda state,t: WaveSymplectic(*state,dt) # evolution rule\n", ")))" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.655221Z", "iopub.status.busy": "2024-04-30T09:15:17.655144Z", "iopub.status.idle": "2024-04-30T09:15:17.656912Z", "shell.execute_reply": "2024-04-30T09:15:17.656684Z" } }, "outputs": [], "source": [ "assert norm_infinity(solution-solution2)<1e-12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.4 Using the Hamiltonian class\n", "\n", "We provide a class devoted to Hamiltonian functions, which are common in mathematics. The wave Hamiltonian has a separable structure, with one generic positive quadratic part an an Euclidean part. It can be specified as follows." ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.658237Z", "iopub.status.busy": "2024-04-30T09:15:17.658163Z", "iopub.status.idle": "2024-04-30T09:15:17.696696Z", "shell.execute_reply": "2024-04-30T09:15:17.696291Z" } }, "outputs": [], "source": [ "reload_packages()" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.698359Z", "iopub.status.busy": "2024-04-30T09:15:17.698270Z", "iopub.status.idle": "2024-04-30T09:15:17.699984Z", "shell.execute_reply": "2024-04-30T09:15:17.699746Z" } }, "outputs": [], "source": [ "Wave = QuadraticHamiltonian(Q_mat,1) # 1 is for the identity matrix" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.701303Z", "iopub.status.busy": "2024-04-30T09:15:17.701220Z", "iopub.status.idle": "2024-04-30T09:15:17.703240Z", "shell.execute_reply": "2024-04-30T09:15:17.703009Z" } }, "outputs": [], "source": [ "assert np.abs(Wave.H(u_cont,u_disc) - WaveHamiltonian(u_cont,u_disc,c_positive,dx)) < 1e-11" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we solve Hamilton's ODE with the symplectic Euler scheme, in the $q$ variant where the position is updated first." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.704607Z", "iopub.status.busy": "2024-04-30T09:15:17.704530Z", "iopub.status.idle": "2024-04-30T09:15:17.710103Z", "shell.execute_reply": "2024-04-30T09:15:17.709857Z" } }, "outputs": [], "source": [ "solution3 = Wave.integrate(u_cont,0.*u_cont,scheme='Euler-q',niter=200,T=200*dt,path=True)\n", "solution3 = np.stack(solution3[:2],axis=1).T\n", "assert norm_infinity(solution-solution3) < 1e-12" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Other schemes are available as well, such as the Verlet second order symplectic scheme." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.711443Z", "iopub.status.busy": "2024-04-30T09:15:17.711364Z", "iopub.status.idle": "2024-04-30T09:15:17.717523Z", "shell.execute_reply": "2024-04-30T09:15:17.717307Z" } }, "outputs": [], "source": [ "solV = Wave.integrate(u_cont,0.*u_smooth,scheme='Verlet-p',niter=200,T=200*dt,path=True)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:17.718907Z", "iopub.status.busy": "2024-04-30T09:15:17.718824Z", "iopub.status.idle": "2024-04-30T09:15:19.907637Z", "shell.execute_reply": "2024-04-30T09:15:19.907317Z" } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "animation_curve(X,solV[0].T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we note that the Hamiltonian class can take care of assembling the sparse matrices, if needed and if the hamiltonian is separable and known to be quadratic." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:19.910689Z", "iopub.status.busy": "2024-04-30T09:15:19.910566Z", "iopub.status.idle": "2024-04-30T09:15:19.914381Z", "shell.execute_reply": "2024-04-30T09:15:19.914131Z" } }, "outputs": [], "source": [ "Hq = lambda q : Q(q,q,c_positive,dx)/2.\n", "Hp = lambda p : I(p,p)/2.\n", "Wave2 = QuadraticHamiltonian(Hq,Hp)\n", "Wave2.set_spmat(u_cont) # Argument just needs to be correctly shaped" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "execution": { "iopub.execute_input": "2024-04-30T09:15:19.915868Z", "iopub.status.busy": "2024-04-30T09:15:19.915784Z", "iopub.status.idle": "2024-04-30T09:15:19.917603Z", "shell.execute_reply": "2024-04-30T09:15:19.917357Z" } }, "outputs": [], "source": [ "assert np.abs(Wave.H(u_cont,u_disc) - Wave2.H(u_cont,u_disc)) < 1e-11" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Format de la Cellule Texte Brut", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.8" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }