{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "133e372e-98c8-43d7-9aa2-fef97b714f5c",
      "metadata": {},
      "source": [
        "# Physics-Informed Neural Network (PINN)\n",
        "**Marc BUFFAT**, mechanical department, University Lyon 1\n",
        "\n",
        "\n",
        "![(C) Zhang et al 2026, Nature](images/PINNs.png)\n",
        "\n",
        "*(C) Zhang, P., Zhang, H., Zhou, J. et al., Sci Rep. 16, 12760, 2026 Natural*\n",
        "\n",
        "*Certain content was co-constructed with the support of a generative AI tool then reread and validated by the teacher.*"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b312889a-3473-468f-af6b-61002b276efe",
      "metadata": {},
      "outputs": [],
      "source": [
        "from validation.libIA_GPU import Init_torchGPU\n",
        "\n",
        "try: cuda_dev\n",
        "except NameError: cuda_dev = Init_torchGPU(4, 0)\n",
        "if cuda_dev is None: cuda_dev = \"cpu\" \n",
        "print(cuda_dev)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "27186d94-009d-4e18-abf1-64fd16cdb86e",
      "metadata": {},
      "source": [
        "**Physics-Informed Neural Networks (PINN)** are a method which consists of training a neural network not only on data, but also on **physical laws** described by differential equations. The idea was popularized by Maziar Raissi and his collaborators in 2019 {cite:ts}`raissi2019`.\n",
        "\n",
        "**references**\n",
        "- http://neuralnetworksanddeeplearning.com/chap4.html\n",
        "- https://karpathy.github.io/2019/04/25/recipe/"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f13f3ea5-3b9e-487e-b722-c742f92362f1",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "\n",
        "print(torch.cuda.is_available())\n",
        "device = torch.device(cuda_dev)\n",
        "print(\"CUDA device:\",cuda_dev)\n",
        "print(\"Torch CUDA device:\",device,\" threads:\",torch.get_num_threads(),torch.get_num_interop_threads())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d7921d49-42ef-4ed1-9cda-99e26864d460",
      "metadata": {},
      "source": [
        "## Function approximation\n",
        "Function approximation is one of the fundamental applications of neural networks. With PyTorch, it's simple to build a network that can learn a function from data.\n",
        "\n",
        "Let's take an example where we want to approximate the function:\n",
        "$f(x)=sin⁡(x)$\n",
        "\n",
        "We define the data `xdata` and `ydata`\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "adee5802-2cbc-4b64-bca6-a6a41a32b008",
      "metadata": {},
      "outputs": [],
      "source": [
        "import torch\n",
        "import torch.nn as nn\n",
        "import matplotlib.pyplot as plt\n",
        "# Training data: attention without gradient\n",
        "Nd = 100\n",
        "xdata = torch.linspace(-2 * torch.pi, 2 * torch.pi, Nd,device=device).reshape(-1, 1)\n",
        "ydata = torch.sin(xdata)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "abb70c37-feef-42fe-9b1f-f0d57595fc6b",
      "metadata": {},
      "outputs": [],
      "source": [
        "# neural network\n",
        "class NeuralNet(nn.Module):\n",
        "    def __init__(self):\n",
        "        super().__init__()\n",
        "        self.model = nn.Sequential(\n",
        "            nn.Linear(1, 32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32, 32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32, 1)\n",
        "        )\n",
        "    def forward(self, x):\n",
        "        return self.model(x)\n",
        "\n",
        "model = NeuralNet().to(device)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "36ad3f3f-36f3-4f63-933a-7de15365b8f5",
      "metadata": {},
      "outputs": [],
      "source": [
        "# cost function\n",
        "criterion = nn.MSELoss()\n",
        "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "131d7fea-e50b-4a1a-9dac-86a1f7d3a38d",
      "metadata": {},
      "outputs": [],
      "source": [
        "# network training\n",
        "num_trainable = sum(p.numel() for p in model.parameters() if p.requires_grad)\n",
        "print(f\"Model Trainable parameters: {num_trainable}\")\n",
        "epochs = 3000\n",
        "for epoch in range(epochs):\n",
        "    optimizer.zero_grad()\n",
        "    ypred = model(xdata)\n",
        "    loss = torch.mean((ypred-ydata)**2)\n",
        "    loss.backward()\n",
        "    optimizer.step()\n",
        "    if epoch % 200 == 0:\n",
        "        print(f\"Epoch {epoch} : Loss = {loss.item():.6f}\")\n",
        "print(f\"Erreur finale = {loss.item():.6f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a2f5ded5-2445-4ea8-8de0-2cb7c29cf49f",
      "metadata": {},
      "source": [
        "### Validation\n",
        "To validate, we define collocation points `x` (note: `xdata` and `x` are not necessarily identical).\n",
        "\n",
        "If we want to calculate the derivatives of the function $f(x)$, we define the collocation points x with the option `x.requires_grad=True`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0407b48e-7223-4b48-b59b-82fad80b1bc0",
      "metadata": {},
      "outputs": [],
      "source": [
        "# collocation points where we define the function\n",
        "N = 2*Nd\n",
        "x = torch.linspace(-2 * torch.pi, 2 * torch.pi, N,device=device).reshape(-1, 1)\n",
        "# if we want to calculate the derivatives\n",
        "x.requires_grad=True\n",
        "# test the model\n",
        "with torch.no_grad():\n",
        "    y_pred = model(x)\n",
        "xn = xdata.cpu()\n",
        "yn = ydata.cpu()\n",
        "xp = x.detach().cpu()\n",
        "yp = y_pred.cpu()\n",
        "plt.figure(figsize=(8,8))\n",
        "plt.subplot(2,1,1)\n",
        "plt.plot(xn, yn, label=\"sin(x)\")\n",
        "plt.plot(xp, yp, '--', label=\"Approximation NN\")\n",
        "plt.legend()\n",
        "plt.grid()\n",
        "plt.subplot(2,1,2)\n",
        "yp = model(xdata).detach().cpu()\n",
        "plt.plot(xn, yn - yp)\n",
        "plt.title(\"Erreur\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eb84dbb6-dc00-4df5-99c4-09ec25104a57",
      "metadata": {},
      "source": [
        "### Calculation of the derivative\n",
        "If we calculated the approximation with `x.requires_grad=True`, we can calculate the derivative with the function `torch.autograd.grad`, which uses automatic differentiation in PyTorch."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "462a6034-d6f5-4714-ac15-eb2614653be4",
      "metadata": {},
      "outputs": [],
      "source": [
        "ypred = model(x)\n",
        "dypred = torch.autograd.grad(\n",
        "        ypred,\n",
        "        x,\n",
        "        grad_outputs=torch.ones_like(ypred),\n",
        "        create_graph=True\n",
        "    )[0]\n",
        "dyp = dypred.detach().cpu()\n",
        "dye = np.cos(xp)\n",
        "# trace\n",
        "plt.figure(figsize=(8,8))\n",
        "plt.subplot(2,1,1)\n",
        "plt.title(\"Dérivée de la fonction $f(x)=sin(x)$\")\n",
        "plt.plot(xp,dyp,label='dérivée PINN')\n",
        "plt.plot(xp,dye,label='exacte')\n",
        "plt.legend()\n",
        "plt.grid()\n",
        "plt.subplot(2,1,2)\n",
        "plt.plot(xp, dye - dyp)\n",
        "plt.title(\"Erreur\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ceb3cf15-5b3d-486e-a28b-38293022ce0b",
      "metadata": {},
      "source": [
        "### universal approximation theorem\n",
        "\n",
        "The universal approximation theorem [Wikipedia](https://fr.wikipedia.org/wiki/Th%C3%A9or%C3%A8me_d%27approximation_universelle) states that a neural network with a sufficiently large hidden layer (with a nonlinear activation function like Tanh, ReLU or Sigmoid) can approximate any continuous function arbitrarily well on a compact domain.\n",
        "\n",
        "**Principle of nonlinear combinations**\n",
        "\n",
        "  - Each neuron shifts, stretches or contracts the input in a slightly different way, then passes it through a non-linear activation function to produce a partial response.\n",
        "  - These partial responses are strongly active in some regions and almost inactive in others, like small fragments of function.\n",
        "  - The output layer then combines these fragments by weighted sum to form the overall function.\n",
        "  - Ultimately, the neural network does not describe a general rule in a single block: it constructs a complex function by stacking many small responses.\n",
        "\n",
        "**Meaning of approximation**\n",
        "\n",
        "  - The universal approximation theorem does not mean that the network always reproduces the original function exactly.\n",
        "  - Its precise meaning is as follows: if an error tolerance is fixed in advance, it is possible to reduce the gap between functions below this tolerance.\n",
        "  - The central idea is therefore not **exactly identical**, but **as close as desired**.\n",
        "  - This distinction makes it possible to avoid confusing representation capacity and real precision.\n",
        "\n",
        "$$ \\hat{f}(x)= \\sum_{i=1}^N w_i \\sigma(a_i x_i + b_i) $$\n",
        "with\n",
        "$$ \\max_{x\\in K} |f(x)-\\hat{f}(x)| < \\epsilon$$\n",
        "\n",
        "**importance and limitation**\n",
        "\n",
        "The universal approximation theorem provides a solid mathematical basis showing that a neural network is not only a simple linear model, but can also represent extremely complex continuous functions. It is therefore often cited as a theoretical starting point to explain why Deep Learning can address a wide variety of problems. \n",
        "\n",
        "However, this theorem only guarantees that the representation is possible: it does not mean that this representation will be effective with few neurons, nor that learning will be simple and fast. In real problems, the number of required neurons can become very large, and difficulties related to optimization, data quantity and generalization arise separately. In other words, **this theorem is essential for understanding the potential of neural networks, but it does not automatically guarantee their performance in practice.**\n",
        "\n",
        "Excerpt from the ZeroMathAI program.\n",
        "- https://zeromathai.com/fr/universal-approximation-theorem-fr/"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "06e4a83a-ee54-4440-93b2-357e0b14879c",
      "metadata": {},
      "source": [
        "### Comparison of approximation\n",
        "\n",
        "We can compare the neural network approximation with a cubic spline approximation.\n",
        "\n",
        "In the case of the sine function, an approximation with 20 splines gives a better result than with the neural network with 1153 parameters!\n",
        "\n",
        "**It is therefore not necessary to systematically use this neural network approximation**."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b0063c8a-bbce-4424-8123-b0fa4dd04cc1",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from scipy.interpolate import CubicSpline\n",
        "\n",
        "N = 20\n",
        "x = np.linspace(-2 * np.pi, 2 * np.pi, N)\n",
        "y = np.sin(x)\n",
        "cs = CubicSpline(x, y)\n",
        "xs = torch.linspace(-2 * np.pi, 2 * np.pi,200)\n",
        "ys = np.sin(xs)\n",
        "ypred = cs(xs)\n",
        "plt.figure(figsize=(8,8))\n",
        "plt.subplot(2,1,1)\n",
        "plt.title(\"Approximation classique\")\n",
        "plt.plot(xs,ys, label=\"sin(x)\")\n",
        "plt.plot(xs,ypred, '--', label=\"Splines\")\n",
        "plt.legend()\n",
        "plt.grid()\n",
        "plt.subplot(2,1,2)\n",
        "plt.plot(xs, ys - ypred)\n",
        "plt.title(\"Erreur\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "369b9197-27f8-43bf-94e4-e1ddd1aac510",
      "metadata": {},
      "source": [
        "**Practical advice**\n",
        "\n",
        "- Normalize the data if the inputs or outputs have very different scales.\n",
        "- Reserve a validation set to monitor overfitting.\n",
        "- Choose the network capacity (number of layers and neurons) based on the complexity of the function to be approximated.\n",
        "- Use torch.utils.data.DataLoader for large datasets.\n",
        "- Monitor the loss level on the training and validation sets in order to adjust the hyper parameters (learning rate, network size, number of epochs).\n",
        "\n",
        "This diagram forms the basis of many regression applications with PyTorch, whether to approximate an analytical function, a physical model or an empirical relationship from experimental data."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "aacc56c2-b1e0-42f5-a64c-7b27dadee7fc",
      "metadata": {},
      "source": [
        "## General principle of a PINN to solve an ODE\n",
        "\n",
        "Suppose we want to solve a differential equation\n",
        "\n",
        "$$\n",
        "\\mathcal{N}[u(x)] = 0,\n",
        "$$\n",
        "\n",
        "with boundary conditions or initial conditions.\n",
        "\n",
        "Instead of discretizing the domain (as with finite difference or finite element methods), we directly approximate the solution using a neural network:\n",
        "\n",
        "$$\n",
        "u(x) \\approx u_\\theta(x),\n",
        "$$\n",
        "\n",
        "where $\\theta$ represents the network weights.\n",
        "\n",
        "The objective is then to find the parameters $\\theta$ which make this function compatible with physical laws.\n",
        "\n",
        "---\n",
        "\n",
        "### Using automatic differentiation\n",
        "\n",
        "One of the key points of PINNs is that the derivatives are calculated using the **automatic differentiation** of PyTorch or TensorFlow.\n",
        "\n",
        "For example, if $u_\\theta = u_\\theta(x),$\n",
        "then $\\frac{\\partial u_\\theta}{\\partial x}$\n",
        "is obtained automatically by\n",
        "\n",
        "```python\n",
        "du_dx = torch.autograd.grad(\n",
        "    outputs=u,\n",
        "    inputs=x,\n",
        "    grad_outputs=torch.ones_like(u),\n",
        "    create_graph=True\n",
        ")[0]\n",
        "```\n",
        "\n",
        "This calculation is not based on a numerical approximation of the derivatives, but on an **automatic numerical differentiation** of the approximation by neural network.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "81cceb2e-f0d7-43a5-b505-21af93a50ca3",
      "metadata": {},
      "source": [
        "### Calculation of the physical residue\n",
        "\n",
        "Consider the ODE\n",
        "\n",
        "$$\n",
        "\\frac{du}{dt}+u=0.\n",
        "$$\n",
        "\n",
        "The neural network predicts $u_\\theta(t)$. We then calculate the derivative:\n",
        "\n",
        "$$\n",
        "\\frac{du_\\theta}{dt},\n",
        "$$\n",
        "\n",
        "then the **residual** of the equation\n",
        "\n",
        "$$\n",
        "r(t)=\\frac{du_\\theta}{dt}+u_\\theta.\n",
        "$$\n",
        "\n",
        "If $u_\\theta(t)$ is solution of the ODE, then $r(t)=0 \\;\\forall t$\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f9508d59-4222-4f7f-b1c4-1ed2d21babf5",
      "metadata": {},
      "source": [
        "### Calculation Cost Function\n",
        "\n",
        "The training consists of minimizing several terms:\n",
        "\n",
        "$$ \\mathcal{L} = \n",
        "\\mathcal L_{\\mathrm{phys}}\n",
        "+\n",
        "\\mathcal L_{\\mathrm{BC}}\n",
        "+\n",
        "\\mathcal L_{\\mathrm{data}}.\n",
        "$$\n",
        "\n",
        "**1. Physical loss**\n",
        "\n",
        "It requires the function $u_\\theta(t)$ to verify the differential equation at certain points.\n",
        "\n",
        "$$\n",
        "\\mathcal L_{\\mathrm{phys}} =\n",
        "\\frac1N\n",
        "\\sum_{i=1}^{N}\n",
        "r(x_i)^2.\n",
        "$$\n",
        "\n",
        "The $x_i$ are called **collocation points**.\n",
        "\n",
        "---\n",
        "\n",
        "**2. Initial or boundary conditions**\n",
        "\n",
        "For example to impose the initial condition $u(0)=1.$, we add\n",
        "\n",
        "$$\n",
        "\\mathcal L_{\\mathrm{BC}} = \n",
        "(u_\\theta(0)-1)^2.\n",
        "$$\n",
        "\n",
        "---\n",
        "\n",
        "**3. Experimental data (optional)**\n",
        "\n",
        "If we have some measurements\n",
        "$(x_i,u_i),$\n",
        "we add\n",
        "\n",
        "$$\n",
        "\\mathcal L_{\\mathrm{data}} = \n",
        "\\frac1M\n",
        "\\sum\n",
        "(u_\\theta(x_i)-u_i)^2.\n",
        "$$\n",
        "\n",
        "Thus, PINNs can work:\n",
        "\n",
        "* without data (physics alone guides learning);\n",
        "* with some data (hybrid case);\n",
        "* with a lot of data (reinforcing physical coherence).\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b857bb43-d109-44f3-b459-5cb9d1948dd7",
      "metadata": {},
      "source": [
        "### Principle of the algorithm\n",
        "\n",
        "```text\n",
        "            Choix des points x\n",
        "                  │\n",
        "                  ▼\n",
        "         Réseau de neurones\n",
        "              uθ(x)\n",
        "                  │\n",
        "        ┌─────────┴──────────┐\n",
        "        ▼                    ▼\n",
        " Conditions          Dérivées automatiques\n",
        "initiales                 (autograd)\n",
        "        │                    │\n",
        "        └─────────┬──────────┘\n",
        "                  ▼\n",
        "          Résidu de l'équation\n",
        "                  │\n",
        "                  ▼\n",
        "        Fonction de coût totale\n",
        "                  │\n",
        "                  ▼\n",
        "     Descente de gradient (Adam/L-BFGS)\n",
        "                  │\n",
        "                  ▼\n",
        "       Mise à jour des paramètres θ\n",
        "```\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "50a25cee-8009-41ed-94db-8d8b4a19e21d",
      "metadata": {},
      "source": [
        "### Why does this work?\n",
        "\n",
        "The **universal approximation theorem** guarantees that a sufficiently large neural network can approximate a very large class of continuous functions. PINNs exploit this property to search for the function that best satisfies the physical equations and the imposed conditions.\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "56f11af3-1e4b-420e-9a72-2a3f24655380",
      "metadata": {},
      "source": [
        "### Advantages\n",
        "\n",
        "* No complex mesh to build.\n",
        "* Derivatives are obtained with high precision thanks to automatic differentiation.\n",
        "* Easy to extend to multidimensional problems.\n",
        "* Naturally integrate physical knowledge.\n",
        "* Can solve **inverse problems**, where unknown parameters are estimated along with the solution.\n",
        "\n",
        "---\n",
        "\n",
        "### Limitations\n",
        "\n",
        "* Training can be slow, especially for large problems.\n",
        "*Equations that are very steep, strongly nonlinear or with discontinuities are difficult to process.\n",
        "* The choice of weights between the different components of the loss function is delicate.\n",
        "* For many industrial PDEs, classic numerical methods (finite elements, finite volumes, finite differences) often remain faster and more precise.\n",
        "\n",
        "**In summary**, a PINN replaces the traditional numerical solver with a neural network trained to directly satisfy physics equations. Solving an ODE or a PDE then becomes an optimization problem where the network parameters are adjusted to minimize the residual of the equation and respect the initial or boundary conditions."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "47527048-a582-42b1-9ffe-756460964708",
      "metadata": {},
      "source": [
        "## Example 1: ODE of order 1\n",
        "Here is a complete example of a **Physics-Informed Neural Network (PINN)** in **PyTorch** to solve an ordinary differential equation (ODE).\n",
        "\n",
        "Consider the following ODE:\n",
        "\n",
        "$$\n",
        "\\frac{dy}{dt} = -y, \\qquad y(0)=1\n",
        "$$\n",
        "\n",
        "The analytical solution is:\n",
        "\n",
        "$$\n",
        "y(t)=e^{-t}.\n",
        "$$\n",
        "\n",
        "The principle of a PINN is to simultaneously minimize:\n",
        "\n",
        "* the residue of the differential equation,\n",
        "* the initial (or boundary) conditions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2ba0b8d1-75a5-47af-8be5-452a7bd71e43",
      "metadata": {},
      "outputs": [],
      "source": [
        "import torch\n",
        "import torch.nn as nn\n",
        "import matplotlib.pyplot as plt\n",
        "# Device choice\n",
        "print(\"use of\",device)\n",
        "# \n",
        "# Neural network\n",
        "# \n",
        "class PINN(nn.Module):\n",
        "    def __init__(self):\n",
        "        super().__init__()\n",
        "        self.net = nn.Sequential(\n",
        "            nn.Linear(1, 32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32, 32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32, 32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32, 1)\n",
        "        )\n",
        "\n",
        "    def forward(self, x):\n",
        "        return self.net(x)\n",
        "\n",
        "model = PINN().to(device)\n",
        "# \n",
        "# Collocation points\n",
        "# \n",
        "N = 100\n",
        "t = torch.linspace(0, 5, N).view(-1,1).to(device)\n",
        "t.requires_grad = True\n",
        "# Initial condition\n",
        "t0 = torch.tensor([[0.0]], device=device)\n",
        "y0 = torch.tensor([[1.0]], device=device)\n",
        "optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)\n",
        "# \n",
        "# Training\n",
        "# \n",
        "epochs = 5000\n",
        "total_params = sum(p.numel() for p in model.parameters() if p.requires_grad)\n",
        "print(f'Total number of parameters: {total_params}')\n",
        "\n",
        "for epoch in range(epochs):\n",
        "    optimizer.zero_grad()\n",
        "    # Prediction\n",
        "    y = model(t)\n",
        "    # Derivative dy/dt\n",
        "    dy_dt = torch.autograd.grad(\n",
        "        outputs=y,\n",
        "        inputs=t,\n",
        "        grad_outputs=torch.ones_like(y),\n",
        "        create_graph=True\n",
        "    )[0]\n",
        "    # ODE residue\n",
        "    residual = dy_dt + y\n",
        "    loss_ode = torch.mean(residual**2)\n",
        "    # Initial condition\n",
        "    pred0 = model(t0)\n",
        "    loss_ic = (pred0 - y0).pow(2).mean()\n",
        "    # Cost function\n",
        "    loss = loss_ode + loss_ic\n",
        "    # descent\n",
        "    loss.backward()\n",
        "    optimizer.step()\n",
        "    # residue display\n",
        "    if epoch % 500 == 0:\n",
        "        print(f\"Epoch {epoch:5d}  Loss={loss.item():.3e}\")\n",
        "\n",
        "# \n",
        "# Evaluation and solution plotting\n",
        "# \n",
        "with torch.no_grad():\n",
        "    t_test = torch.linspace(0,5,200).view(-1,1).to(device)\n",
        "    y_pred = model(t_test).cpu()\n",
        "\n",
        "t_exact = t_test.cpu()\n",
        "y_exact = torch.exp(-t_exact)\n",
        "\n",
        "plt.figure(figsize=(8,5))\n",
        "plt.plot(t_exact, y_exact, label=\"Solution exacte\")\n",
        "plt.plot(t_exact, y_pred, \"--\", label=\"PINN\")\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"y\")\n",
        "plt.legend()\n",
        "plt.grid()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8aacf797-0fb1-4dcc-9426-13cdf21b81f9",
      "metadata": {},
      "source": [
        "### Explanation\n",
        "\n",
        "The PINN learns a function $y_\\theta(t)$ parameterized by a neural network.\n",
        "\n",
        "At each iteration:\n",
        "\n",
        "1. The predicted network (y(t)).\n",
        "2. PyTorch automatically calculates the derivative\n",
        "   $\\frac{dy_\\theta}{dt}$\n",
        "   thanks to `torch.autograd.grad`.\n",
        "3. We construct the residue of the ODE:\n",
        "   $r(t)=\\frac{dy_\\theta}{dt}+y_\\theta.$\n",
        "4. The loss function is\n",
        "   $$\n",
        "   \\mathcal{L} =\n",
        "   \\underbrace{\\frac1N\\sum_i r(t_i)^2}*{\\text{equation}}\n",
        "   +\n",
        "   \\underbrace{\\left(y*\\theta(0)-1\\right)^2}_{\\text{initial condition}}.\n",
        "   $$\n",
        "\n",
        "The network is thus trained to simultaneously satisfy the differential equation and the initial condition."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7a653f76-51df-491a-84fe-5cdab02618ac",
      "metadata": {},
      "source": [
        "### Generalization\n",
        "\n",
        "For any ODE\n",
        "\n",
        "$$\n",
        "\\frac{dy}{dt}=f(t,y),\n",
        "$$\n",
        "\n",
        "just replace\n",
        "\n",
        "```python\n",
        "residual = dy_dt + y\n",
        "```\n",
        "\n",
        "by\n",
        "\n",
        "```python\n",
        "residual = dy_dt - f(t, y)\n",
        "```\n",
        "\n",
        "where `f(t, y)` is a Python function defining the right-hand side."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ece8fd1f-b841-43d3-a800-e8e354c3c8a2",
      "metadata": {},
      "source": [
        "## Example 2: the harmonic oscillator\n",
        "\n",
        "$$\n",
        "\\frac{d^2y}{dt^2}+y=0,\n",
        "$$\n",
        "\n",
        "with initial conditions\n",
        "\n",
        "$$\n",
        "y(0)=1,\\qquad y'(0)=0.\n",
        "$$\n",
        "\n",
        "The exact solution is\n",
        "\n",
        "$$\n",
        "y(t)=\\cos(t).\n",
        "$$\n",
        "\n",
        "The idea is to have this solution learned by a PINN.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f168b4fc-0613-4c30-a1d3-564ef7dd8441",
      "metadata": {},
      "outputs": [],
      "source": [
        "import torch\n",
        "import torch.nn as nn\n",
        "import matplotlib.pyplot as plt\n",
        "# Device choice\n",
        "print(\"using\",device)\n",
        "# =========================================================\n",
        "# Neural network\n",
        "# =========================================================\n",
        "class PINN(nn.Module):\n",
        "    def __init__(self):\n",
        "        super().__init__()\n",
        "\n",
        "        self.net = nn.Sequential(\n",
        "            nn.Linear(1,64),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(64,64),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(64,64),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(64,1)\n",
        "        )\n",
        "\n",
        "    def forward(self,x):\n",
        "        return self.net(x)\n",
        "\n",
        "model = PINN().to(device)\n",
        "# =========================================================\n",
        "# Collocation points\n",
        "# =========================================================\n",
        "N = 200\n",
        "t = torch.linspace(0,10,N).reshape(-1,1).to(device)\n",
        "t.requires_grad = True\n",
        "\n",
        "optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)\n",
        "# =========================================================\n",
        "# Training\n",
        "# =========================================================\n",
        "epochs = 5000\n",
        "total_params = sum(p.numel() for p in model.parameters() if p.requires_grad)\n",
        "print(f'Total number of parameters: {total_params}')\n",
        "\n",
        "for epoch in range(epochs):\n",
        "    optimizer.zero_grad()\n",
        "    # prediction\n",
        "    y = model(t)\n",
        "    # first derivative\n",
        "    dy = torch.autograd.grad(\n",
        "        y,\n",
        "        t,\n",
        "        grad_outputs=torch.ones_like(y),\n",
        "        create_graph=True\n",
        "    )[0]\n",
        "    # second derivative\n",
        "    d2y = torch.autograd.grad(\n",
        "        dy,\n",
        "        t,\n",
        "        grad_outputs=torch.ones_like(dy),\n",
        "        create_graph=True\n",
        "    )[0]\n",
        "    # ODE residue\n",
        "    residual = d2y + y\n",
        "    loss_phys = torch.mean(residual**2)\n",
        "    # =========================================================\n",
        "    # Initial conditions\n",
        "    # =========================================================\n",
        "    t0 = torch.tensor([[0.0]], requires_grad=True,device=device)\n",
        "    y0 = model(t0)\n",
        "    dy0 = torch.autograd.grad(\n",
        "        y0,\n",
        "        t0,\n",
        "        grad_outputs=torch.ones_like(y0),\n",
        "        create_graph=True\n",
        "    )[0]\n",
        "    loss_ic = (y0-1)**2 + dy0**2\n",
        "    # =========================================================\n",
        "    loss = loss_phys + loss_ic\n",
        "    loss.backward()\n",
        "    optimizer.step()\n",
        "    if epoch % 500 == 0:\n",
        "        print(epoch, loss.item())\n",
        "\n",
        "# =========================================================\n",
        "# verification\n",
        "# =========================================================\n",
        "with torch.no_grad():\n",
        "    t_test = torch.linspace(0,10,300,device=device).view(-1,1)\n",
        "    y_pred = model(t_test).cpu()\n",
        "t_exact = t_test.cpu()\n",
        "y_exact = np.cos(t_exact)\n",
        "\n",
        "plt.figure(figsize=(8,5))\n",
        "plt.plot(t_exact,y_exact,label=\"Exact\")\n",
        "plt.plot(t_exact,y_pred,\"--\",label=\"PINN\")\n",
        "plt.legend()\n",
        "plt.grid()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f7996812-f03c-4eb4-9f86-12ce335d65da",
      "metadata": {},
      "source": [
        "### Explanation\n",
        "\n",
        "For a second order ODE, the network directly approximates the function:\n",
        "\n",
        "$$\n",
        "y(t)\\approx y_\\theta(t).\n",
        "$$\n",
        "\n",
        "The derivatives are calculated by automatic differentiation:\n",
        "\n",
        "$$\n",
        "\\frac{dy_\\theta}{dt},\n",
        "\\qquad\n",
        "\\frac{d^2y_\\theta}{dt^2}.\n",
        "$$\n",
        "\n",
        "In PyTorch, this is done with two successive calls to `torch.autograd.grad`:\n",
        "\n",
        "```python\n",
        "dy = torch.autograd.grad(\n",
        "    y, t,\n",
        "    grad_outputs=torch.ones_like(y),\n",
        "    create_graph=True\n",
        ")[0]\n",
        "\n",
        "d2y = torch.autograd.grad(\n",
        "    dy, t,\n",
        "    grad_outputs=torch.ones_like(dy),\n",
        "    create_graph=True\n",
        ")[0]\n",
        "```\n",
        "\n",
        "The residue of the equation is then\n",
        "\n",
        "$$\n",
        "r(t)=\\frac{d^2y_\\theta}{dt^2}+y_\\theta,\n",
        "$$\n",
        "\n",
        "and the physical loss is\n",
        "\n",
        "$$\n",
        "L_{\\mathrm{phys}} =\n",
        "\\frac1N\\sum_i r(t_i)^2.\n",
        "$$\n",
        "\n",
        "The initial conditions are added in the form of penalties:\n",
        "\n",
        "$$\n",
        "L_{\\mathrm{IC}}=\n",
        "(y_\\theta(0)-1)^2\n",
        "+\n",
        "\\left(\\frac{dy_\\theta}{dt}(0)\\right)^2.\n",
        "$$\n",
        "\n",
        "The total cost function is\n",
        "\n",
        "$$\n",
        "L = L_{\\mathrm{phys}} + L_{\\mathrm{IC}}.\n",
        "$$\n",
        "\n",
        "### Generalization\n",
        "\n",
        "For a second order ODE of the form\n",
        "\n",
        "$$\n",
        "y'' = f(t,y,y'),\n",
        "$$\n",
        "\n",
        "just replace the residue with:\n",
        "\n",
        "```python\n",
        "residual = d2y - f(t, y, dy)\n",
        "```\n",
        "\n",
        "where `f` is a Python function. This structure is valid for most second-order ODEs, whether linear or nonlinear. It can then be extended to ODE systems by extracting several components from the network and constructing a residual for each equation."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ca54c70b-aea0-4d0a-a3c7-41dee47f42a8",
      "metadata": {},
      "source": [
        "## Example 3: Reverse Problem\n",
        "\n",
        "**Inverse problems** are one of the most interesting applications of PINNs. Unlike the direct problem (where the parameters of the equation are known), the goal here is to **estimate unknown physical parameters** from a few measurements.\n",
        "\n",
        "Consider the following ODE:\n",
        "\n",
        "$$\n",
        "\\frac{dy}{dt} = -\\lambda y,\n",
        "$$\n",
        "\n",
        "with $y(0)=1,$ where **$\\lambda$** is unknown.\n",
        "\n",
        "Suppose that we have some noisy measurements of $y(t)$ and that we wish to find $\\lambda$.\n",
        "\n",
        "---\n",
        "\n",
        "### Principle\n",
        "\n",
        "The network approximates the solution\n",
        "$y_\\theta(t).$\n",
        "\n",
        "But this time, the physical parameter is also optimized:\n",
        "\n",
        "```python\n",
        "lambda_param = nn.Parameter(torch.tensor([0.5]))\n",
        "```\n",
        "\n",
        "It is added to the list of trainable parameters."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "df183b80-2b66-43a6-b5f7-b6072ac1817e",
      "metadata": {},
      "outputs": [],
      "source": [
        "import torch\n",
        "import torch.nn as nn\n",
        "import matplotlib.pyplot as plt\n",
        "# Device choice\n",
        "print(\"using\",device)\n",
        "# =========================================================\n",
        "# Data generation\n",
        "# =========================================================\n",
        "lambda_true = 2.0\n",
        "t_data = torch.linspace(0, 2, 20,device=device).reshape(-1,1)\n",
        "y_data = torch.exp(-lambda_true*t_data)\n",
        "# noise\n",
        "y_data += 0.01*torch.randn_like(y_data)\n",
        "# =========================================================\n",
        "# PINN\n",
        "# =========================================================\n",
        "class PINN(nn.Module):\n",
        "\n",
        "    def __init__(self):\n",
        "        super().__init__()\n",
        "        self.net = nn.Sequential(\n",
        "            nn.Linear(1,32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32,32),\n",
        "            nn.Tanh(),\n",
        "            nn.Linear(32,1)\n",
        "        )\n",
        "\n",
        "    def forward(self,x):\n",
        "        return self.net(x)\n",
        "\n",
        "model = PINN().to(device)\n",
        "# =========================================================\n",
        "# Unknown parameter\n",
        "# =========================================================\n",
        "lambda_param = nn.Parameter(torch.tensor([0.5],device=device))\n",
        "# =========================================================\n",
        "# Optimizer\n",
        "# =========================================================\n",
        "optimizer = torch.optim.Adam(\n",
        "    list(model.parameters())+[lambda_param],\n",
        "    lr=1e-3\n",
        ")\n",
        "# =========================================================\n",
        "# Collocation points\n",
        "# =========================================================\n",
        "t_phys = torch.linspace(0,2,100,device=device).reshape(-1,1)\n",
        "t_phys.requires_grad=True\n",
        "\n",
        "# =========================================================\n",
        "# Training\n",
        "# =========================================================\n",
        "total_params = sum(p.numel() for p in model.parameters() if p.requires_grad)\n",
        "print(f'Total number of parameters: {total_params}')\n",
        "\n",
        "for epoch in range(5000):\n",
        "    optimizer.zero_grad()\n",
        "    # \n",
        "    # Physical\n",
        "    # \n",
        "    y = model(t_phys)\n",
        "    dy = torch.autograd.grad(\n",
        "        y,\n",
        "        t_phys,\n",
        "        grad_outputs=torch.ones_like(y),\n",
        "        create_graph=True\n",
        "    )[0]\n",
        "    residual = dy + lambda_param*y\n",
        "    loss_phys = torch.mean(residual**2)\n",
        "    # \n",
        "    # Data\n",
        "    # \n",
        "    y_pred = model(t_data)\n",
        "    loss_data = torch.mean((y_pred-y_data)**2)\n",
        "    # \n",
        "    # Initial condition\n",
        "    # \n",
        "    t0 = torch.tensor([[0.0]],device=device)\n",
        "    loss_ic = (model(t0)-1.)**2\n",
        "    # \n",
        "    loss = loss_phys + loss_data + loss_ic\n",
        "    loss.backward()\n",
        "    optimizer.step()\n",
        "    if epoch%500==0:\n",
        "        print(epoch,\n",
        "              loss.item(),\n",
        "              lambda_param.item())\n",
        "\n",
        "print(\"\\nEstimated value:\", lambda_param.item(), \"exacte:\",lambda_true)\n",
        "\n",
        "# =========================================================\n",
        "# Comparison\n",
        "# =========================================================\n",
        "with torch.no_grad():\n",
        "    t_test = torch.linspace(0,2,200,device=device).reshape(-1,1)\n",
        "    y_pred = model(t_test).cpu()\n",
        "\n",
        "t = t_test.cpu()\n",
        "plt.plot(t,np.exp(-lambda_true*t),label=\"Exacte\")\n",
        "plt.plot(t,y_pred,'--',label=\"PINN\")\n",
        "plt.scatter(t_data.cpu(),y_data.cpu(),color='red')\n",
        "plt.legend()\n",
        "plt.grid()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2de00832-ae46-4e88-976a-7e06aa469984",
      "metadata": {},
      "source": [
        "### Explanation\n",
        "\n",
        "The network simultaneously learns:\n",
        "\n",
        "* the function $y(t)$,\n",
        "* the $\\lambda$ parameter.\n",
        "\n",
        "The loss function is written\n",
        "\n",
        "$$\n",
        "L=\n",
        "L_{\\text{phys}}\n",
        "+\n",
        "L_{\\text{data}}\n",
        "+\n",
        "L_{\\text{IC}}\n",
        "$$\n",
        "\n",
        "with\n",
        "\n",
        "1. Physical Loss\n",
        "\n",
        "$$\n",
        "L_{\\text{phys}}\n",
        "=\n",
        "\\frac1N\n",
        "\\sum_i\n",
        "\\left(\n",
        "\\frac{dy_\\theta}{dt}\n",
        "+\\lambda y_\\theta\n",
        "\\right)^2.\n",
        "$$\n",
        "\n",
        "2. Loss of Data\n",
        "\n",
        "$$\n",
        "L_{\\text{data}}=\n",
        "\\frac1M\n",
        "\\sum_i\n",
        "\\left(\n",
        "y_\\theta(t_i)-y_i\n",
        "\\right)^2.\n",
        "$$\n",
        "\n",
        "3. Loss on the Initial Condition\n",
        "\n",
        "$$\n",
        "L_{\\text{IC}}=\n",
        "(y_\\theta(0)-1)^2.\n",
        "$$\n",
        "\n",
        "---\n",
        "\n",
        "### Why does this work?\n",
        "\n",
        "The loss gradient depends on **two types of parameters**:\n",
        "\n",
        "* the weights of the neural network $\\theta$;\n",
        "* the physical parameter $\\lambda$.\n",
        "\n",
        "As $\\lambda$ is declared as a `nn.Parameter`, PyTorch automatically calculates\n",
        "\n",
        "$$\n",
        "\\frac{\\partial L}{\\partial \\lambda},\n",
        "$$\n",
        "\n",
        "and the optimizer updates $\\lambda$ exactly like the network weights.\n",
        "\n",
        "---\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "67ac413a-0ab0-4767-9eb2-69ef0acf5353",
      "metadata": {},
      "source": [
        "## Real-world applications of PINNs\n",
        "\n",
        "This approach is used to identify physical parameters from limited data, for example:\n",
        "\n",
        "* diffusion coefficient in a heat equation;\n",
        "* viscosity of a fluid in the Navier–Stokes equations;\n",
        "* Young's modulus or thermal conductivity of a material;\n",
        "* reaction constants in chemical kinetics;\n",
        "* parameters of an epidemiological model (for example the transmission rate in an SIR model);\n",
        "* damping and stiffness parameters in mechanical systems.\n",
        "\n",
        "It is precisely in this context where there are few measurements but a reliable physical model that PINNs often offer their greatest value.\n",
        "\n",
        "### Possible extensions\n",
        "\n",
        "This structure easily extends to more complex problems:\n",
        "\n",
        "* **ODE systems** (Lotka-Volterra, pendulum, Lorenz, etc.) by bringing several components out of the network.\n",
        "* **Higher order differential equations**, calculating successive derivatives with `torch.autograd.grad`.\n",
        "* **Partial differential equations (PDE)** such as the heat, Burgers or Navier-Stokes equation, using multiple input variables (`x`, `t`, etc.) and corresponding partial derivatives.\n",
        "* **Inverse PINNs**, where some physical parameters of the ODE are also learned by the network by declaring them as optimizable parameters (`torch.nn.Parameter`)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "46774e34-1cd0-4363-b577-fc3ecc5e1e7d",
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "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.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}