Second order ODE
Problem setup
We will solve an ODE:
with initial conditions
For \(t \in [0,0.25]\). The exact solution is \(y(t)=\frac{50}{81}+\frac{5}{9}t + \frac{31}{81} e^{9t} - 2e^t\).
Implementation
This description goes through the implementation of a solver for the above ODE step-by-step. First, the DeepXDE module is imported:
import deepxde as dde
import numpy as np
We begin by defining a computational geometry. We can use a built-in class TimeDomain as follows
geom = dde.geometry.TimeDomain(0, 0.25)
Next, we express the residual of the ODE:
def ode(t, y):
dy_dt = dde.grad.jacobian(y, t)
d2y_dt2 = dde.grad.hessian(y, t)
return d2y_dt2 - 10 * dy_dt + 9 * y - 5 * t
The first argument to ode is the network input, i.e., the \(t\)-coordinate. The second argument is the network output, i.e., the solution \(y(t)\), but here we use y as the name of the variable.
We define the initial condition, setting the value of the function at \(t=0\) to -1.
ic1 = dde.icbc.IC(geom, lambda x: -1, lambda _, on_initial: on_initial)
Now we deal with the initial condition \(y'(0)=2\).
The location of the intial condition is defined by a simple Python function. The function should return True for those points satisfying \(t=0\) and False otherwise (note that because of rounding-off errors, it is often wise to use dde.utils.isclose to test whether two floating point values are equivalent). In this function, the argument t to boundary is the network input and is a \(d\)-dim vector, where \(d\) is the dimension and \(d=1\) in this case. Then a boolean on_boundary is used as the second argument. If the point t (the first argument) is on the boundary of the geometry, then on_boundary is True, otherwise, on_boundary is False.
def boundary_l(t, on_boundary):
return on_boundary and dde.utils.isclose(t[0], 0)
Now we define a function that returns the error of the initial condition, \(y'(0)=2\), which is the difference between the derivative of the output of the network at 0, and 2. The function takes arguments (inputs, outputs, X) and outputs a tensor of size N x 1, where N is the length of inputs. inputs and outputs are the network input and output tensors, respectively; X is the NumPy array of the inputs.
def error_2(inputs, outputs, X):
return dde.grad.jacobian(outputs, inputs, i=0, j=None) - 2
Then, the initial condition is defined by
ic2 = dde.icbc.OperatorBC(geom, error_2, boundary_l)
Now, we have specified the geometry, PDE residual, and the initial conditions. We then define the PDE problem. Note: If you use X in func, then do not set num_test when you define dde.data.PDE or dde.data.TimePDE, otherwise DeepXDE would throw an error. In this case, the training points will be used for testing, and this will not affect the network training and training loss. This is a bug of DeepXDE, which cannot be fixed in an easy way for all backends.
data = dde.data.TimePDE(geom, ode, [ic1, ic2], 16, 2, solution=func, num_test=500)
The number 16 is the number of training residual points sampled inside the domain, and the number 2 is the number of training points sampled on the boundary. The argument solution=func is the reference solution to compute the error of our solution, and can be ignored if we don’t have a reference solution. We use 500 residual points for testing the PDE residual.
Next, we choose the network. Here, we use a fully connected neural network of depth 4 (i.e., 3 hidden layers) and width 50:
layer_size = [1] + [50] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)
Now, we have the PDE problem and the network. We build a Model, choose the optimizer, set the learning rate to 0.001, and train the network for 15000 iterations. We set the weight of the ODE loss to 0.01, and the weights of the two ICs to 1. We also compute the \(L^2\) relative error as a metric during training.
model = dde.Model(data, net)
model.compile("adam", lr=.001, loss_weights=[0.01, 1, 1], metrics=["l2 relative error"])
losshistory, train_state = model.train(iterations=15000)
Complete code
"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch, paddle"""
import deepxde as dde
import numpy as np
def ode(t, y):
dy_dt = dde.grad.jacobian(y, t)
d2y_dt2 = dde.grad.hessian(y, t)
return d2y_dt2 - 10 * dy_dt + 9 * y - 5 * t
def func(t):
return 50 / 81 + t * 5 / 9 - 2 * np.exp(t) + (31 / 81) * np.exp(9 * t)
geom = dde.geometry.TimeDomain(0, 0.25)
def boundary_l(t, on_initial):
return on_initial and dde.utils.isclose(t[0], 0)
def bc_func1(inputs, outputs, X):
return outputs + 1
def bc_func2(inputs, outputs, X):
return dde.grad.jacobian(outputs, inputs, i=0, j=None) - 2
ic1 = dde.icbc.IC(geom, lambda x: -1, lambda _, on_initial: on_initial)
ic2 = dde.icbc.OperatorBC(geom, bc_func2, boundary_l)
data = dde.data.TimePDE(geom, ode, [ic1, ic2], 16, 2, solution=func, num_test=500)
layer_size = [1] + [50] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile(
"adam", lr=0.001, metrics=["l2 relative error"], loss_weights=[0.01, 1, 1]
)
losshistory, train_state = model.train(iterations=10000)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)