TP: Physics Informed Neural Networks (PINNs)¶
I) Equation différentielle ordinaire¶
Dans ce TP nous commencerons par voir comment simuler numériquement une EDO à l'aide d'un PINN. On considère l'EDO suivante. $$ y'(t) + a \, y(t) = f(t), $$ avec $y(0) = 0$, et $$ f(t) = \begin{cases} e^{-bt}\sin(t) & \text{si } t \geq 0, \\ 0 & \text{si } t < 0, \end{cases} $$ qui admet pour solution explicite: $$ y(t) = \begin{cases} \frac{1}{(a-b)^2 + 1} \big( ((a-b) \sin(t) - \cos(t)) e^{-bt} + e^{-at} \big) & \text{si } t \geq 0, \\ 0 & \text{si } t < 0. \end{cases} $$
In [3]:
import numpy as np
import matplotlib.pyplot as plt
In [4]:
# paramètres de l'EDO
a = 0.1
b = 0.1
# grille de points auxquelles on veut évaluer la solution y
t = np.linspace(-5., 30., 100)
# définition de la solution explicite
sol_expl = lambda t: ( ( (a-b) * np.sin(t) - np.cos(t) ) * np.exp(-b*t) + np.exp(-a*t) ) / ((a-b)**2 + 1) *(t>0)
In [5]:
plt.plot(t, sol_expl(t), label='solution y(t)')
plt.plot(t, np.sin(t)*np.exp(-b*t)*(t>0) , label='source f(t)')
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.title('explicite solution and source f', fontsize=18)
plt.show()
Résolution de l'EDO avec le solver odeint¶
Pour résoudre des EDO avec python, on peut utiliser la fonction odeint de la librairie scipy.integrate. Pour cela on définit la fonction suivante.
In [6]:
from scipy.integrate import odeint
In [7]:
def eq(u, t, a, b):
dudt = -a * u + np.sin(t)*np.exp(-b*t)*(t>0)
return dudt
In [8]:
# grille de points sur laquelle le solver odeint évalue la solution
t_eval = np.linspace(-5, 30, 50)
# condition initiale de l'EDO
y0 = 0.
# calcul de la solution
sol = odeint(eq, y0, t_eval, args=(a,b))
le calcul de la solution est quasi instantanée. Ca ne sera pas le cas avec un réseau de neurones comme on le verra plus loin.
In [9]:
plt.plot(t, sol_expl(t), label='explicite solution')
plt.plot(t_eval, sol, label='odeint sol')
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('explicite solution vs odeint simulation', fontsize=18)
plt.show()
Sans surprise l'approximation est excellente.
Approche classique d'apprentisage¶
- Commençons par supposer qu'on a des observations de la solution de l'EDO, et approchons ces observations à l'aide d'un réseau de neurones.
- Une fois le réseau entrainé, nous verrons :
- comment celui-ci généralise
- et s'il vérifie le modèle mathématique ou pas.
In [18]:
# nombre d'observation
n_obs = 70
# temps auxquelles on fait les observations
t_obs = np.linspace(t[0], t[-1], n_obs)
# observation de la solutions
y_obs = sol_expl(t_obs) + 0.0*np.random.randn(n_obs) # on peut éventuellement rajouter du bruit
In [9]:
# on rajoute la condition initiale qui est une donnée du problème
t_obs = np.concatenate(([0.],t_obs))
y_obs = np.concatenate(([y0],y_obs))
In [10]:
plt.plot(t, sol_expl(t), label='ground truth')
plt.plot(t_obs, y_obs, 'o', label='data')
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('observations of the ode', fontsize=18)
plt.show()
Question: Entrainer un réseau de neurone qui, à partir des données, permet d'approcher les observations de la solution de l'équation différentielle. Illustrer numériquement les valeurs des résidus¶
$$ \mathcal{N}'_\theta(t) + a \, \mathcal{N}_\theta(t) - f(t),\qquad t\in[-5,30]$$
où $\mathcal{N}_\theta$ est le réseau de neurones entrainé sur les données. Que peut-on observer?¶
Définissons une classe de réseau de neurones.
In [11]:
import torch
import torch.nn as nn
import torch.optim as optim
In [11]:
class MLP(nn.Module):
def __init__(self, layers, activations):
"""
Parameters:
- layers: List[int], the number of neurons in each layer (including input and output layers).
- activations: List[nn.Module], the activation functions for each layer.
The length of activations should be len(layers) - 2 (one for each hidden layer).
"""
super(MLP, self).__init__()
if len(activations) != len(layers) - 2:
raise ValueError("The number of activation functions must match the number of hidden layers.")
# Build the network
self.layers = nn.Sequential()
for i in range(len(layers) - 1):
self.layers.add_module(f"linear_{i}", nn.Linear(layers[i], layers[i + 1]))
if i < len(activations): # Apply activation to all but the last layer
self.layers.add_module(f"activation_{i}", activations[i])
def forward(self, x):
return self.layers(x)
In [13]:
# Liste contenant les caractéristiques des couches (profondeur et nombre de neurones
# Input, hidden layers, output
layers = [1, 50, 50, 50, 1]
# Liste contenant les fonctions d'activations pour les couches cachées
activations = [nn.Tanh(), nn.Tanh(), nn.Tanh()] # Different activations for each hidden layer
In [14]:
#définition du modèle
model = MLP(layers, activations)
print(model)
MLP(
(layers): Sequential(
(linear_0): Linear(in_features=1, out_features=50, bias=True)
(activation_0): Tanh()
(linear_1): Linear(in_features=50, out_features=50, bias=True)
(activation_1): Tanh()
(linear_2): Linear(in_features=50, out_features=50, bias=True)
(activation_2): Tanh()
(linear_3): Linear(in_features=50, out_features=1, bias=True)
)
)
In [15]:
# fonction d'entrainement
def train(model, optimizer, loss_fct, t_data, y_data, epochs=1000):
"""
Trains a given model on specific data for a defined number of epochs.
Parameters:
- model : nn.Module
The neural network model to be trained.
- optimizer : torch.optim.Optimizer
The optimizer used to update the model's weights.
- loss_fct :
The loss function to measure the error between the model's predictions and the target data.
- t_data : torch.Tensor
Input data (time variables).
- y_data : torch.Tensor
Target data corresponding to the input data.
- epochs : int, optional (default 1000)
The number of epochs (iterations over the entire dataset) for training.
"""
for epoch in range(epochs):
# reset gradients for the optimizer
optimizer.zero_grad()
# compare the model to the data
y_pred = model(t_data)
loss = loss_fct(y_pred, y_data)
# optimisation step
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print(f"Epoch {epoch}, Loss: {loss.item()}")
print(f"Epoch {epoch}, Loss: {loss.item()}")
In [16]:
# définition de la méthode d'optimisation
optimizer = optim.Adam(model.parameters(), lr=0.01)
In [17]:
# fonction de coût quadratique
loss_fct = torch.nn.MSELoss(reduction='sum')
In [18]:
# mise en forme des données en torch.tensor
t_data = torch.tensor(t_obs, dtype=torch.float32).reshape(-1,1)
y_data = torch.tensor(y_obs, dtype=torch.float32).reshape(-1,1)
In [19]:
# entrainement
train(model, optimizer, loss_fct, t_data, y_data, 30000)
Epoch 0, Loss: 12.267515182495117 Epoch 1000, Loss: 0.34444481134414673 Epoch 2000, Loss: 0.11400552093982697 Epoch 3000, Loss: 0.01016367319971323 Epoch 4000, Loss: 0.017028529196977615 Epoch 5000, Loss: 0.00033501762663945556 Epoch 6000, Loss: 0.0002265641960548237 Epoch 7000, Loss: 0.06964896619319916 Epoch 8000, Loss: 0.014804917387664318 Epoch 9000, Loss: 0.1380312293767929 Epoch 10000, Loss: 0.009359326213598251 Epoch 11000, Loss: 0.005478056147694588 Epoch 12000, Loss: 0.01858140155673027 Epoch 13000, Loss: 0.0006765498546883464 Epoch 14000, Loss: 0.02274276688694954 Epoch 15000, Loss: 0.008032940328121185 Epoch 16000, Loss: 0.010999591089785099 Epoch 17000, Loss: 0.06218697130680084 Epoch 18000, Loss: 0.011863981373608112 Epoch 19000, Loss: 0.0015520744491368532 Epoch 20000, Loss: 0.00029264503973536193 Epoch 21000, Loss: 0.00029000442009419203 Epoch 22000, Loss: 0.028675494715571404 Epoch 23000, Loss: 0.0026121537666767836 Epoch 24000, Loss: 0.2914021611213684 Epoch 25000, Loss: 0.030146511271595955 Epoch 26000, Loss: 0.00021879335690755397 Epoch 27000, Loss: 0.001125902053900063 Epoch 28000, Loss: 0.02337130531668663 Epoch 29000, Loss: 5.272788985166699e-05 Epoch 29999, Loss: 0.0010860570473596454
In [20]:
# évaluation du model sur une grille de temps fine
t_val = torch.tensor(t, dtype=torch.float32).reshape(-1,1)
# comparaison de la solution explicite avec les prédictions du réseau de neurones
plt.plot(t, sol_expl(t), label="ground truth")
plt.plot(t_obs, y_obs, 'o', label="data")
plt.plot(t_val, model(t_val).detach(), label="predictions") # on utilise .detach()"
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('fit of the ode data', fontsize=18)
plt.show()
L'adéquation de la solution explicite et du réseau sont très bonne même en des points qui n'ont pas servi à l'entrainement. Cependant, le reseau de neurone vérifie-t-il l'EDO, respecte-t-il la physique?
In [21]:
t_val = t_val.requires_grad_()
y = model(t_val)
In [22]:
# calcul des dérivées en temps du réseau sur la grille
y_t = torch.autograd.grad(y, t_val, grad_outputs=torch.ones_like(y))[0]
In [23]:
# évaluation de l'équation mathématique pour notre réseau
residual = y_t + a * y - torch.sin(t_val)*torch.exp(-b*t_val) * (t_val>0)
torch.sum(residual**2)
Out[23]:
tensor(6.7147, grad_fn=<SumBackward0>)
Si le réseau vérifie l'EDO, residual ne doit avoir que des petites valeurs.
In [24]:
with torch.no_grad():
plt.plot(t_val, residual)
plt.xlabel('t', fontsize=18)
plt.ylabel("y'(t)+ay(t)-f(t)", fontsize=18)
plt.title('residual', fontsize=18)
plt.show()
C'est raté! Notre réseau ne vérifie pas l'EDO. Ces prédictions non donc pas de valeur par rapport au modèle mathématique.
Si on augmente la taille des données d'entrainement les residus
$$ \mathcal{N}'_\theta(t) + a \, \mathcal{N}_\theta(t) - f(t)$$
sont plus petits. Avec beaucoup de données, un réseau de neurones peut apprendre l'équation.
Essayons maintenant d'entrainer notre modèle avec une loss contenant les informations sur l'EDO.
Physics Informed Neural Networks (PINNs)¶
On va entrainer un réseau de neurones avec une loss faisant intervenir les termes de l'EDO: \begin{align*} L(\theta) &= L_{edo}(\theta) + L_{ci}(\theta) \\ & = \frac{1}{N} \sum_{j=1}^N |\mathcal{N}'_\theta(t_j) + a \mathcal{N}_\theta(t_j) - f(t_j)|^2 + (\mathcal{N}_\theta(t = 0) - y_0)^2. \end{align*}
On reprend exactement la même architecture que précédemment.
In [15]:
layers = [1, 50, 50, 50, 1] # Input, hidden layers, output
activations = [nn.Tanh(), nn.Tanh(), nn.Tanh()] # Different activations for each hidden layer
# Initialize the model
y0 = 0.
model_pinn = MLP(layers, activations)
# Print the model architecture
print(model_pinn)
MLP(
(layers): Sequential(
(linear_0): Linear(in_features=1, out_features=50, bias=True)
(activation_0): Tanh()
(linear_1): Linear(in_features=50, out_features=50, bias=True)
(activation_1): Tanh()
(linear_2): Linear(in_features=50, out_features=50, bias=True)
(activation_2): Tanh()
(linear_3): Linear(in_features=50, out_features=1, bias=True)
)
)
Question: Entrainer un PINN permettant de simuler numériquement la solution de l'équation différentielle. Illustrer numériquement les résidus¶
$$ \mathcal{N}'_\theta(t) + a \, \mathcal{N}_\theta(t) - f(t),\qquad t\in[-5,30]$$
où $\mathcal{N}_\theta$ est maintenant un PINN. Que peut-on observer par rapport aux résidus du réseau uniquement entrainé sur les données? On prendra autant de point de collocation que de points d'observations dans le cas d'entrainement sur des données.¶
On définit la fonction d'entrainement du réseau
In [16]:
def train_pinn(model, optimizer, t, y0, a, f, epochs=1000):
for epoch in range(epochs):
"""
Parameters:
- model (torch.nn.Module): The neural network model that approximates the solution to the ODE.
- optimizer (torch.optim.Optimizer): The optimizer used to update the model's weights during training.
- t (torch.Tensor): A tensor representing the input time values for which the ODE solution is being approximated.
- y0 (torch.Tensor): The initial condition for the ODE, specifying the solution value at t = 0.
- a (float): A coefficient in the ODE, representing the linear term multiplier in the equation y'+ay−f=0.
- f (torch.Tensor): The forcing function in the ODE.
- epochs (int, optional): The number of training iterations (default: 1000).
"""
# reset gradients for the optimizer
optimizer.zero_grad()
### compute the loss involving the ODE terms
# evaluate the model on the ODE grid
y = model(t)
# compute the derivatives involved in the ode
y_t = torch.autograd.grad(y, t, grad_outputs=torch.ones_like(y), create_graph=True)[0]
residual = y_t + a * y - f
# evaluate the residual, that is the mismatch with the ode
loss_ode = residual.square().sum()
#evaluate the mismatch with the initial condition
loss_ic = (model(torch.tensor([0.]))-y0)**2
#compute the loss
loss = loss_ode + loss_ic
# optimization step
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print(f"Epoch {epoch}, Loss: {loss.item()}")
print(f"Epoch {epoch}, Loss: {loss.item()}")
On reprend la même grille que pour l'approche ML standard, mais dans l'approche PINNs les points sont appelés points de collocations
In [19]:
# nombre d'observation
n_colloc = 70
# On définit la grille de temps pour L_{ode} sur laquelle on entraine le réseau
t_colloc = np.linspace(t[0], t[-1], n_colloc)
# observation de la solutions
y_colloc = sol_expl(t_colloc) + 0.0*np.random.randn(n_obs) # on peut éventuellement rajouter du bruit
In [20]:
# conversion en torch.tensor
tt_colloc = torch.tensor(t_colloc, dtype=torch.float32).reshape(-1,1).requires_grad_()
# on calcul y' par rapport à ces valeurs
# et le terme source
with torch.no_grad():
f = np.sin(tt_colloc) * np.exp(-b*tt_colloc) *(tt_colloc>0)
In [21]:
# on définit la méthode d'optimisation
optimizer = optim.Adam(model_pinn.parameters(), lr=0.001)
In [22]:
train_pinn(model_pinn, optimizer, tt_colloc, y0, a, f, 20000)
Epoch 0, Loss: 4.799078464508057 Epoch 1000, Loss: 0.8097332715988159 Epoch 2000, Loss: 0.2179691344499588 Epoch 3000, Loss: 0.13749445974826813 Epoch 4000, Loss: 0.06717421859502792 Epoch 5000, Loss: 0.05131564661860466 Epoch 6000, Loss: 0.007739121560007334 Epoch 7000, Loss: 0.00018009135965257883 Epoch 8000, Loss: 0.00010567763820290565 Epoch 9000, Loss: 0.00044300477020442486 Epoch 10000, Loss: 0.000339063408318907 Epoch 11000, Loss: 4.1931303712772205e-05 Epoch 12000, Loss: 3.615720197558403e-05 Epoch 13000, Loss: 4.1350034734932706e-05 Epoch 14000, Loss: 0.00012688837887253612 Epoch 15000, Loss: 3.369563273736276e-05 Epoch 16000, Loss: 0.00017000740626826882 Epoch 17000, Loss: 3.424451278988272e-05 Epoch 18000, Loss: 2.7771031454904005e-05 Epoch 19000, Loss: 6.0767855757148936e-05 Epoch 19999, Loss: 7.017493044259027e-05
In [24]:
# évaluation du model sur une grille de temps fine
t_val = torch.tensor(t, dtype=torch.float32).reshape(-1,1)
# comparaison de la solution explicite avec les prédictions du réseau de neurones
plt.plot(t, sol_expl(t), label="ground truth")
plt.plot(t_colloc, y_colloc, 'o', label="data")
plt.plot(t_val, model_pinn(t_val).detach(), label="predictions") # on utilise .detach()"
plt.legend(fontsize=14)
plt.xlabel('t', fontsize=18)
plt.ylabel('y(t)', fontsize=18)
plt.title('fit of the ode data with a pinn', fontsize=18)
plt.show()
Regardons si le réseaux vérifie l'EDO. Si oui, residual ne doit avoir que des petites valeurs.
In [25]:
t_val = t_val.requires_grad_()
y = model_pinn(t_val)
y_t = torch.autograd.grad(y, t_val, grad_outputs=torch.ones_like(y), create_graph=True)[0]
residual = y_t + a * y - torch.sin(t_val)*torch.exp(-b*t_val) * (t_val>0)
torch.sum(residual**2).item()
Out[25]:
0.0019436399452388287
In [26]:
with torch.no_grad():
plt.plot(t_val, residual)
plt.xlabel('t', fontsize=18)
plt.ylabel("y'(t)+ay(t)-f(t)", fontsize=18)
plt.title('residual', fontsize=18)
plt.show()
C'est beaucoup mieux qu'avec l'approche standard de ML
II) Equation de la chaleur avec condition de Dirichlet sur $[0,1]$¶
On considère l'équation de la chaleur sur $[0,1]$ $$ \frac{\partial u}{\partial t} - \sigma^2\frac{\partial^2 u}{\partial x^2} = 0, \qquad t, x \in [0,1], $$ où
- l'inconnue est une fonction de 2 variables, $u(t,x)$ représentant la température au temps $t$ et au point $x$;
- $\sigma^2$ est le coefficient de diffusion.
À cette équation, on rajoute :
- une condition initiale (température initiale) : $$ u(t=0,x) = f(x) = e^{-40(x-0.5)^2}, \qquad x \in [0,1]; $$
- les conditions de Dirichlet aux bords du domaine : $$ u(t,x=0) = 0, \quad u(t,x=1) = 0, \quad t \geq 0. $$
Les cellules suivantes permettent de calculer une solution approchée à partir d'une décomposition en série de Fourier $$ u(t,x) = \sum_{j=0}^{+\infty} a_j \sin(j \pi x) e^{-(j \pi)^2 \sigma^2 t} \qquad \text{avec} \qquad a_j = 2\int_0^1 f(x) \sin(j \pi x) dx. $$
Remarque: Pour améliorer la simulation numérique de l'équation de la chaleur avec un PINN, on ne met pas les conditions aux bords dans la loss, on considère plutôt l'approximation $$ u(t,x) \simeq x(1-x) \mathcal{N}_\theta(t,x) $$ où $\mathcal{N}_\theta$ est un réseau de neurones à entrainer. La multiplication par la fonction $x\mapsto x(1-x)$ nous garantie que les conditions de Dirichlet seront respectées.
In [166]:
import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt
In [167]:
import torch
import torch.nn as nn
import torch.optim as optim
from torch.autograd import grad
In [168]:
M = 1000
t = np.linspace(0, 1, M)
x = np.linspace(0, 1, M)
sig2 = 0.3
Question: Illustrer numériquement la solution approchée de l'équation de la chaleur par série de Fourier.¶
In [169]:
def f(x):
return np.exp(-40*(x-0.5)**2)
In [170]:
# number of Fourier modes
N = 30
# compute the Fourier coefficients
A = np.zeros(N+1)
for j in range(1,N+1):
A[j], _ = quad(lambda x: np.sin(j*x*np.pi) * f(x) , 0, 1)
A[j] *= 2
In [171]:
# Generate x and t values
x = np.linspace(0, 1, 100)
t = np.linspace(0, 1, 100)
xx, tt = np.meshgrid(x, t)
In [172]:
sol = np.zeros((100, 100))
for j in range(N):
sol += A[j] * np.sin(j*xx * np.pi) * np.exp(-(j*np.pi)**2 * sig2 * tt)
In [173]:
# indices des temps pour dessiner la solution
idx_t = [0, 20, 40, 60, 80, 99]
for idx in idx_t:
plt.plot(x, sol[idx,:], label=f't={t[idx]:.2f}')
plt.xlabel("x", fontsize=18)
plt.ylabel("u(t, x)", fontsize=18)
plt.legend(fontsize=14)
# Adjust layout and show the plot
plt.suptitle(f"Dirichlet boundary conditions ($\sigma^2$ = {sig2})", fontsize=16)
plt.show()
Question: A l'aide d'un PINN simuler numériquement l'équation de la chaleur¶
- On définira une fonction d'entrainement adaptée au problème;
- On commencera par considérer une grille de points de collocation uniforme;
- Ensuite on considérera à chaque epoch une grille en $t$ et $x$ aléatoire de 256 points
- On comparera sur un graphique:
- les prédictions du réseau de neurones après entrainement;
- l'approximation numérique qu'on a calculé à la question précédente.
In [174]:
class PINN_heat(nn.Module):
def __init__(self, layers, activations):
"""
Parameters:
- layers: List[int], the number of neurons in each layer (including input and output layers).
- activations: List[nn.Module], the activation functions for each layer.
The length of activations should be len(layers) - 2 (one for each hidden layer).
"""
super().__init__()
if len(activations) != len(layers) - 2:
raise ValueError("The number of activation functions must match the number of hidden layers.")
self.layers = nn.Sequential()
for i in range(len(layers) - 1):
self.layers.add_module(f"linear_{i}", nn.Linear(layers[i], layers[i + 1]))
if i < len(layers) - 2:
self.layers.add_module(f"activation_{i}", activations[i])
def forward(self, t, x):
inputs = torch.hstack((t, x))
N = self.layers(inputs)
output = x * (1-x) * N # hard code the Dirichlet boundary conditions with the fct x -> x(1-x)
return output
In [175]:
def train_heat(model, optimizer, sig2, batch_size=128, epochs=1000):
for epoch in range(epochs):
t_tmp = torch.rand((batch_size,1)).requires_grad_()
x_tmp = torch.rand((batch_size,1)).requires_grad_()
optimizer.zero_grad()
u = model(t_tmp,x_tmp)
u_t = grad(u, t_tmp, grad_outputs=torch.ones_like(u), create_graph=True)[0]
u_x = grad(u, x_tmp, grad_outputs=torch.ones_like(u), create_graph=True)[0]
u_xx = grad(u_x, x_tmp, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]
residual = u_t - sig2 * u_xx
loss_pde = residual.square().sum()
u0 = model(torch.zeros_like(x_tmp),x_tmp)
loss_ic = (u0 - f_t(x_tmp)).square().sum()
loss = loss_pde + loss_ic
loss.backward()
optimizer.step()
if epoch%1000 == 0:
print(f"Epoch {epoch}, Loss: {loss.item()}")
print(f"Epoch {epoch}, Loss: {loss.item()}")
In [1]:
# condition initial en pytorch
def f_t(x):
return torch.exp(-40*(x-0.5)**2)
In [180]:
learning_rate = 0.001
epochs = 30000
batch_size = 256
In [181]:
layers = [2, 32, 64, 32, 1] # Input, hidden layers, output
activations = [nn.Tanh(), nn.Tanh(), nn.Tanh()] # Different activations for each hidden layer
model_heat_d = PINN_heat(layers, activations)
print(model_heat_d)
PINN_heat(
(layers): Sequential(
(linear_0): Linear(in_features=2, out_features=32, bias=True)
(activation_0): Tanh()
(linear_1): Linear(in_features=32, out_features=64, bias=True)
(activation_1): Tanh()
(linear_2): Linear(in_features=64, out_features=32, bias=True)
(activation_2): Tanh()
(linear_3): Linear(in_features=32, out_features=1, bias=True)
)
)
In [182]:
optimizer_heat_d = torch.optim.Adam(model_heat_d.parameters(), lr=learning_rate)
In [183]:
train_heat(model_heat_d, optimizer_heat_d, sig2, batch_size, epochs)
Epoch 0, Loss: 70.29417419433594 Epoch 1000, Loss: 12.767661094665527 Epoch 2000, Loss: 2.6814143657684326 Epoch 3000, Loss: 1.2109179496765137 Epoch 4000, Loss: 1.4681681394577026 Epoch 5000, Loss: 0.6766965389251709 Epoch 6000, Loss: 0.7916008234024048 Epoch 7000, Loss: 1.274336338043213 Epoch 8000, Loss: 0.29044851660728455 Epoch 9000, Loss: 0.38564619421958923 Epoch 10000, Loss: 1.5259748697280884 Epoch 11000, Loss: 0.1495717465877533 Epoch 12000, Loss: 0.15335559844970703 Epoch 13000, Loss: 0.31482160091400146 Epoch 14000, Loss: 1.4875727891921997 Epoch 15000, Loss: 0.0772721916437149 Epoch 16000, Loss: 0.2222786247730255 Epoch 17000, Loss: 0.12865008413791656 Epoch 18000, Loss: 0.2210649698972702 Epoch 19000, Loss: 0.26046380400657654 Epoch 20000, Loss: 0.3941280245780945 Epoch 21000, Loss: 0.0671190544962883 Epoch 22000, Loss: 0.039129503071308136 Epoch 23000, Loss: 0.15099605917930603 Epoch 24000, Loss: 0.11384232342243195 Epoch 25000, Loss: 0.2939276099205017 Epoch 26000, Loss: 0.13701365888118744 Epoch 27000, Loss: 0.1593925505876541 Epoch 28000, Loss: 0.11049889773130417 Epoch 29000, Loss: 0.09981995820999146 Epoch 29999, Loss: 0.2798042893409729
In [184]:
# on définit les grilles en t et x pour évaluer le réseau
t_plot = torch.tensor(t, dtype=torch.float32)
x_plot = torch.tensor(x, dtype=torch.float32)
tt_plot = torch.tensor(tt, dtype=torch.float32).reshape(-1,1)
xx_plot = torch.tensor(xx, dtype=torch.float32).reshape(-1,1)
In [185]:
xx_plot.shape
Out[185]:
torch.Size([10000, 1])
In [186]:
# Evaluate the model
with torch.no_grad():
u = model_heat_d(tt_plot, xx_plot) #model_heat_d(torch.hstack((t_plot, x_plot)))
In [187]:
# Reshape the output to match the grid
u = u.numpy().reshape(100, 100)
In [188]:
from mpl_toolkits.mplot3d import Axes3D
# Create the 3D plot
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(131, projection='3d')
# Plot the surface
ax.plot_surface(tt, xx, u, cmap='viridis')
# Add labels
ax.set_xlabel('t', fontsize=18)
ax.set_ylabel('x', fontsize=18)
ax.set_zlabel('u(t,x)', fontsize=18)
ax.set_title('PINN Output', fontsize=18)
ax = fig.add_subplot(132, projection='3d')
# Plot the surface
ax.plot_surface(tt, xx, sol, cmap='viridis')
# Add labels
ax.set_xlabel('t', fontsize=18)
ax.set_ylabel('x', fontsize=18)
ax.set_zlabel('u(t,x)', fontsize=18)
ax.set_title('Fourier solution', fontsize=18)
ax = fig.add_subplot(133, projection='3d')
# Plot the surface
ax.plot_surface(tt, xx, u-sol, cmap='viridis')
err = np.max(np.abs(u-sol))/np.max(np.abs(sol)) * 100
# Add labels
ax.set_xlabel('t', fontsize=18)
ax.set_ylabel('x', fontsize=18)
ax.set_title(f'relative error = {err:.2f}%', fontsize=18)
plt.show()
In [189]:
# Define the times to plot the solution
idx_t = [0, 20, 40, 60, 80, 99]
fig, axes = plt.subplots(2, 3, figsize=(15, 8)) # 2 rows, 3 columns for 6 subplots
axes = axes.flatten() # Flatten the axes array for easy indexing
for i, idx in enumerate(idx_t):
axes[i].plot(x, u[idx,:], label='pinn')
axes[i].plot(x, sol[idx,:], label='Fourier')
axes[i].set_title(f"t = {t[idx]:.2f}", fontsize=18)
axes[i].set_xlabel("x", fontsize=18)
axes[i].set_ylabel("u(t, x)", fontsize=18)
axes[i].grid(True)
axes[i].legend(fontsize=10)
# Adjust layout and show the plot
plt.suptitle(f"Heat Equation Solution ($\sigma^2$ = {sig2})", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.95]) # Adjust spacing for the suptitle
plt.show()
