Solving ODEs with SpotOptim’s LinearRegressor
This tutorial demonstrates how to use Physics-Informed Neural Networks (PINNs) to solve ordinary differential equations (ODEs) using SpotOptim’s LinearRegressor class. We’ll solve a first-order ODE with an initial condition, combining data-driven learning with physics constraints.
We want to solve the following ODE:
\[ \frac{dy}{dt} + 0.1 y - \sin\left(\frac{\pi t}{2}\right) = 0 \]
with initial condition:
\[ y(0) = 0 \]
First, let’s import the necessary libraries:
import matplotlib.pyplot as plt
import numpy as np
import torch
import torch.nn as nn
from typing import Tuple
from spotoptim.nn.linear_regressor import LinearRegressor
# Set random seed for reproducibility
torch.manual_seed(123)<torch._C.Generator at 0x112376810>We’ll use SpotOptim’s LinearRegressor class to create a neural network with:
model = LinearRegressor(
input_dim=1,
output_dim=1,
l1=32,
num_hidden_layers=3,
activation="Tanh",
lr=1.0
)
# Get optimizer with custom learning rate
optimizer = model.get_optimizer("Adam", lr=3.0) # 3.0 * 0.001 = 0.003
print(f"Model architecture:")
print(model)Model architecture:
LinearRegressor(
(activation): Tanh()
(network): Sequential(
(0): Linear(in_features=1, out_features=32, bias=True)
(1): Tanh()
(2): Linear(in_features=32, out_features=32, bias=True)
(3): Tanh()
(4): Linear(in_features=32, out_features=32, bias=True)
(5): Tanh()
(6): Linear(in_features=32, out_features=1, bias=True)
)
)We’ll generate exact solution data using a numerical solver (RK2 method):
def oscillator(
n_steps: int = 3000,
t_min: float = 0.0,
t_max: float = 30.0,
y0: float = 0.0,
alpha: float = 0.1,
omega: float = np.pi / 2
) -> Tuple[torch.Tensor, torch.Tensor]:
"""
Solve ODE: dy/dt + alpha*y - sin(omega*t) = 0
using RK2 (midpoint method).
Args:
n_steps: Number of time steps
t_min: Start time
t_max: End time
y0: Initial condition y(t_min)
alpha: Damping coefficient
omega: Forcing frequency
Returns:
Tuple of (time points, solution values) as torch tensors
"""
t_step = (t_max - t_min) / n_steps
t_points = np.arange(t_min, t_min + n_steps * t_step, t_step)[:n_steps]
y = [y0]
for t_current_step_end in t_points[1:]:
t_midpoint = t_current_step_end - t_step / 2.0
y_prev = y[-1]
# Stage 1: intermediate value
slope_at_t_mid = -alpha * y_prev + np.sin(omega * t_midpoint)
y_intermediate = y_prev + (t_step / 2.0) * slope_at_t_mid
# Stage 2: final value
slope_at_t_end = -alpha * y_intermediate + np.sin(omega * t_current_step_end)
y_next = y_prev + t_step * slope_at_t_end
y.append(y_next)
t_tensor = torch.tensor(t_points, dtype=torch.float32).view(-1, 1)
y_tensor = torch.tensor(y, dtype=torch.float32).view(-1, 1)
return t_tensor, y_tensorGenerate the exact solution and sample training data points:
# Generate exact solution (3000 points)
x, y = oscillator()
# Sample training data (every 119th point, giving ~25 points)
x_data = x[0:3000:119]
y_data = y[0:3000:119]
print(f"Exact solution points: {len(x)}")
print(f"Training data points: {len(x_data)}")Exact solution points: 3000
Training data points: 26Visualize the exact solution and training data:
plt.figure(figsize=(10, 6))
plt.plot(x.numpy(), y.numpy(), 'b-', linewidth=2, label="Exact solution y(t)")
plt.scatter(x_data.numpy(), y_data.numpy(), color="tab:orange",
s=50, label="Training data", zorder=5)
plt.xlabel("Time t", fontsize=12)
plt.ylabel("Solution y(t)", fontsize=12)
plt.title("ODE Solution and Training Data", fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Create collocation points where we’ll enforce the physics (ODE) constraints:
# 50 evenly spaced points in [0, 30]
x_physics = torch.linspace(0, 30, 50).view(-1, 1).requires_grad_(True)
print(f"Collocation points: {len(x_physics)}")
print(f"Range: [{x_physics.min().item():.2f}, {x_physics.max().item():.2f}]")Collocation points: 50
Range: [0.00, 30.00]Now we’ll train the neural network using two loss components:
The network predicts values at the training data points:
yh = model(x_data)Mean squared error between predictions and actual data:
loss1 = torch.mean((yh - y_data)**2)The network predicts values at the physics enforcement points:
yhp = model(x_physics)We compute dy/dt using automatic differentiation:
dyhp_dxphysics = torch.autograd.grad(
yhp, x_physics,
torch.ones_like(yhp),
create_graph=True
)[0]The ODE residual at collocation points:
\[ \text{residual} = \frac{dy}{dt} + 0.1 y - \sin\left(\frac{\pi t}{2}\right) \]
physics = dyhp_dxphysics + 0.1 * yhp - torch.sin(np.pi * x_physics / 2)
loss2 = torch.mean(physics**2)Combined loss with weighting factor α:
loss = loss1 + alpha * loss2loss_history_pinn = []
loss2_history_pinn = []
plot_data_points_pinn = []
alpha = 6e-2 # Weight for physics loss
n_epochs = 48000
print("Training Physics-Informed Neural Network...")
print(f"Total epochs: {n_epochs}")
print(f"Physics loss weight (alpha): {alpha}")
for i in range(n_epochs):
optimizer.zero_grad()
# Data Loss: Fit the training data
yh = model(x_data)
loss1 = torch.mean((yh - y_data)**2)
# Physics Loss: Satisfy the ODE at collocation points
yhp = model(x_physics)
dyhp_dxphysics = torch.autograd.grad(
yhp, x_physics,
torch.ones_like(yhp),
create_graph=True
)[0]
physics = dyhp_dxphysics + 0.1 * yhp - torch.sin(np.pi * x_physics / 2)
loss2 = torch.mean(physics**2)
# Total Loss
loss = loss1 + alpha * loss2
loss.backward()
optimizer.step()
# Store history every 100 steps
if (i + 1) % 100 == 0:
loss_history_pinn.append(loss.detach())
loss2_history_pinn.append(loss2.detach())
# Store snapshots for visualization every 10000 steps
if (i + 1) % 10000 == 0:
current_yh_full = model(x).detach()
plot_data_points_pinn.append({
'yh': current_yh_full,
'step': i + 1
})
print(f"Epoch {i+1}/{n_epochs}: Total Loss = {loss.item():.6f}, "
f"Physics Loss = {loss2.item():.6f}")
print("Training completed!")Training Physics-Informed Neural Network...
Total epochs: 48000
Physics loss weight (alpha): 0.06
Epoch 10000/48000: Total Loss = 0.000119, Physics Loss = 0.001705
Epoch 20000/48000: Total Loss = 0.000035, Physics Loss = 0.000473
Epoch 30000/48000: Total Loss = 0.000014, Physics Loss = 0.000225
Epoch 40000/48000: Total Loss = 0.000181, Physics Loss = 0.000391
Training completed!fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Plot total loss
ax1.plot(loss_history_pinn, 'b-', linewidth=1.5)
ax1.set_xlabel('Iteration (×100)', fontsize=11)
ax1.set_ylabel('Total Loss', fontsize=11)
ax1.set_title('Training Progress: Total Loss', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.set_yscale('log')
# Plot physics loss
ax2.plot(loss2_history_pinn, 'r-', linewidth=1.5)
ax2.set_xlabel('Iteration (×100)', fontsize=11)
ax2.set_ylabel('Physics Loss', fontsize=11)
ax2.set_title('Training Progress: Physics Loss (ODE Residual)', fontsize=12)
ax2.grid(True, alpha=0.3)
ax2.set_yscale('log')
plt.tight_layout()
plt.show()
Visualize how the neural network solution evolves during training:
xp_plot = x_physics.detach()
for plot_info in plot_data_points_pinn:
plt.figure(figsize=(10, 6))
# Plot exact solution
plt.plot(x.numpy(), y.numpy(), 'b-', linewidth=2,
label='Exact solution', alpha=0.7)
# Plot training data
plt.scatter(x_data.numpy(), y_data.numpy(), color='tab:orange',
s=50, label='Training data', zorder=5)
# Plot neural network prediction
plt.plot(x.numpy(), plot_info['yh'].numpy(), 'r--', linewidth=2,
label='PINN prediction', alpha=0.8)
# Plot collocation points
plt.scatter(xp_plot.numpy(),
model(xp_plot).detach().numpy(),
color='green', marker='x', s=30,
label='Collocation points', alpha=0.6, zorder=4)
plt.xlabel('Time t', fontsize=12)
plt.ylabel('Solution y(t)', fontsize=12)
plt.title(f'PINN Solution at Epoch {plot_info["step"]}', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Final prediction
model.eval()
with torch.no_grad():
y_final = model(x)
plt.figure(figsize=(12, 7))
# Plot exact solution
plt.plot(x.numpy(), y.numpy(), 'b-', linewidth=2.5,
label='Exact solution', alpha=0.8)
# Plot PINN prediction
plt.plot(x.numpy(), y_final.numpy(), 'r--', linewidth=2,
label='PINN prediction', alpha=0.8)
# Plot training data
plt.scatter(x_data.numpy(), y_data.numpy(), color='tab:orange',
s=80, label='Training data', zorder=5, edgecolors='black', linewidth=0.5)
# Plot collocation points
plt.scatter(x_physics.detach().numpy(),
model(x_physics).detach().numpy(),
color='green', marker='x', s=50,
label='Collocation points', alpha=0.7, zorder=4)
plt.xlabel('Time t', fontsize=13)
plt.ylabel('Solution y(t)', fontsize=13)
plt.title('Final PINN Solution vs Exact Solution', fontsize=15, fontweight='bold')
plt.legend(fontsize=12, loc='best')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Compute absolute error
error = torch.abs(y_final - y)
plt.figure(figsize=(10, 6))
plt.plot(x.numpy(), error.numpy(), 'r-', linewidth=2)
plt.xlabel('Time t', fontsize=12)
plt.ylabel('Absolute Error |y_exact - y_PINN|', fontsize=12)
plt.title('PINN Approximation Error', fontsize=14)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print("\nError Statistics:")
print(f"Maximum absolute error: {error.max().item():.6f}")
print(f"Mean absolute error: {error.mean().item():.6f}")
print(f"Root mean squared error: {torch.sqrt(torch.mean(error**2)).item():.6f}")
Error Statistics:
Maximum absolute error: 0.072458
Mean absolute error: 0.008300
Root mean squared error: 0.013641This tutorial demonstrated how to use SpotOptim’s LinearRegressor class to implement a Physics-Informed Neural Network (PINN) for solving ordinary differential equations. Key takeaways:
The PINN successfully learned to approximate the ODE solution using only sparse training data (~25 points) by incorporating the underlying physics through the differential equation constraint.
The LinearRegressor class integrates seamlessly with SpotOptim for hyperparameter tuning:
from spotoptim import SpotOptim
def train_pinn(X):
"""Objective function for hyperparameter optimization."""
results = []
for params in X:
l1 = int(params[0]) # Hidden layer size
num_layers = int(params[1]) # Number of layers
lr_unified = 10 ** params[2] # Learning rate (log scale)
# Create model
model = LinearRegressor(
input_dim=1, output_dim=1,
l1=l1, num_hidden_layers=num_layers,
activation="Tanh", lr=lr_unified
)
# Train PINN and compute validation error
# ... training code ...
results.append(validation_error)
return np.array(results)
# Optimize hyperparameters
optimizer = SpotOptim(
fun=train_pinn,
bounds=[(16, 128), (1, 5), (-4, 0)],
var_type=["int", "int", "float"],
var_name=["layer_size", "num_layers", "log_lr"],
max_iter=50
)
result = optimizer.optimize()This approach allows you to systematically find the best network architecture and learning rate for your specific PINN problem.