ОДУ второго порядка

Общая информация

Покажем, что ничего страшного в ОДУ второго и более высоких порядков нет - алгоритм всё тот же:

1. Составить выражение для пробной функции, удовлетворяющее заданным условиям вне зависимости от выхода нейронной сети

2. Нужное количество раз взять производную от этого выражения

3. Засунуть всё это в нейронную сеть, предварительно задав её конфигурацию, и долго ждать, пока всё досчитается

Рассмотрим на примере ОДУ второго порядка на отрезке \([a, b]\)

\[ \begin{equation*} \frac{d^{2} f(x)}{d x^{2}}=G\left(x, f(x), \frac{d f(x)}{d x}\right) \end{equation*} \]

Для ОДУ второго порядка уже есть свобода в выборе вида граничных условий.

1) Задача Коши

\[f(a)=A \in \mathbb{R}, \; f^{\prime}(a)=A^{\prime} \in \mathbb{R}\]

Составим пробную функцию в виде

\[ \begin{equation*} \hat{f}(x ; \theta)=A+A^{\prime}(x-a)+(x-a)^{2} N(x ; \theta) \end{equation*} \]

Эта функция удовлетворяет условиям независимо от \(N(x ; \theta)\)

Можно вычислить производные

\[\begin{split} \begin{align*} \frac{\partial \hat{f}(x ; \theta)}{\partial x} & =A^{\prime}+(x-a)\left[2 N(x ; \theta)+(x-a) \frac{\partial N(x ; \theta)}{\partial x}\right] \\ \frac{\partial^{2} \hat{f}(x ; \theta)}{d x^{2}} & =2 N(x ; \theta)+4(x-a) \frac{\partial N(x ; \theta)}{d x}+(x-a)^{2} \frac{\partial^{2} N(x ; \theta)}{d x^{2}} \end{align*} \end{split}\]

\[ \begin{equation*} J(\theta)=\frac{1}{M} \sum_{j=1}^{M}\left(\frac{\partial^{2} \hat{f}\left(x_{j} ; \theta\right)}{\partial x^{2}}-G\left(x_{j}, \hat{f}\left(x_{j} ; \theta\right), \frac{\partial \hat{f}\left(x_{j} ; \theta\right)}{\partial x}\right)\right)^{2} \end{equation*} \]

2) Краевая задача (условия Дирихле)

\[f(a)=A \in \mathbb{R}, \; f(b)=B \in \mathbb{R}\]

\[ \begin{equation*} \hat{f}(x ; \theta)=A \frac{(b-x)}{b-a}+B \frac{(x-a)}{b-a}+(x-a)(b-x) N(x ; \theta) \end{equation*} \]

Остальные увлекательные шаги для условий Дирихле в виде вычисления производных читатель может проделать самостоятельно.

Пример (задача Коши)

Приведём уравнение, точное решение которого будет совпадать с решением примера 2 из пункта про ОДУ первого порядка.

\[ \begin{equation*} \frac{d^{2} f(x)}{d x^{2}}=-\frac{1}{5} \frac{d f(x)}{d x}-f(x)-\frac{1}{5} e^{-\frac{x}{5}} \cos (x) \end{equation*} \]

\[ \begin{equation*} f(0)=0, \; f^{\prime}(0)=1 \end{equation*} \]

Аналитическое решение

\[ \begin{equation*} f(x)=e^{-\frac{x}{5}} \sin (x) \end{equation*} \]

С учётом граничных условий пробная функция запишется как

\[ \begin{equation*} \hat{f}(x ; \theta)=x+x^{2} N(x ; \theta) \end{equation*} \]

import torch
import torch.utils.data
import numpy as np
import matplotlib.pyplot as plt
from torch.autograd import grad
import time

class DataSet(torch.utils.data.Dataset):
    def __init__(self, numSamples, xRange):
        self.dataIn  = torch.linspace(xRange[0], xRange[1], numSamples, requires_grad=True).view(-1,1)
    def __len__(self):
        return len(self.dataIn)

    def __getitem__(self, idx):
        return self.dataIn[idx]

class Fitter(torch.nn.Module):
    def __init__(self, numHiddenNodes):
        super(Fitter, self).__init__()
        self.fc1 = torch.nn.Linear(in_features = 1, out_features = numHiddenNodes)
        self.fc2 = torch.nn.Linear(in_features = numHiddenNodes, out_features = 1)

    def forward(self, x):
        h = torch.tanh(self.fc1(x))
        y = self.fc2(h)
        return y

def train(network, loader, lossFn, optimiser, numEpochs):
    cost_list=[]
    network.train(True)
    for epoch in range(numEpochs):
        for batch in loader:
            n_out = network(batch)
            dndx = grad(n_out, batch, torch.ones_like(n_out), retain_graph=True, create_graph=True)[0]
            d2ndx2 =grad(dndx, batch, torch.ones_like(dndx), retain_graph=True)[0]
            f_trial = trialFunc(batch, n_out)
            df_trial = dTrialFunc(batch, n_out, dndx)
            d2f_trial = d2TrialFunc(batch,n_out,dndx,d2ndx2)
            diff_eq = diffEq(batch, f_trial, df_trial, d2f_trial)
            cost = lossFn(diff_eq, torch.zeros_like(diff_eq))
            cost.backward()
            optimiser.step()
            optimiser.zero_grad()

        cost_list.append(cost.detach().numpy())
    network.train(False)
    return cost_list


def plotNetwork(network, epoch):
    x    = torch.linspace(0, 10, 120, requires_grad=True).view(-1,1)
    N    = network.forward(x)
    f_trial = trialFunc(x, N)
    dndx = grad(N, x, torch.ones_like(N), retain_graph=True, create_graph=True)[0]
    d2ndx2 =grad(dndx, x, torch.ones_like(dndx), retain_graph=True)[0]
    df_trial = dTrialFunc(x, N, dndx)
    d2f_trial = d2TrialFunc(x,N,dndx,d2ndx2)
    diff_eq = diffEq(x, f_trial, df_trial, d2f_trial)
    cost = lossFn(diff_eq, torch.zeros_like(diff_eq))

    exact = solution(x)
    MSECost = lossFn(f_trial, exact)
    print("MSE между аналитическим и полученным решениями = ", MSECost.item())
    exact = exact.detach().numpy()
    x = x.detach().numpy()
    N = N.detach().numpy()
    plt.plot(x, trialFunc(x,N), 'r-', label = "Neural Network Output")
    plt.plot(x, exact, 'b.', label = "True Solution")
    plt.xlabel("x", fontsize = 16)
    plt.ylabel("f(x)", fontsize = 16)
    plt.legend(loc = "upper right", fontsize = 16)
    plt.title(str(epoch) + " Epochs", fontsize = 16)
    plt.show()
    return


def trialFunc(x, n_out):
    return x + (x**2 * n_out)

def dTrialFunc(x, n_out, dndx):
    return 1 + (2*x*n_out) + (x**2 * dndx)

def d2TrialFunc(x,n_out,dndx,d2ndx2):
    return 2*n_out + (4*x*dndx) + (x**2 * d2ndx2)

def diffEq(x, f_trial, df_trial, d2f_trial):
    LHS = d2f_trial + (1/5)*df_trial + f_trial
    RHS = -(1/5) * torch.exp(-x/5) * torch.cos(x)
    return LHS - RHS

def solution(x):
    return torch.exp(-x/5) * torch.sin(x)

network = Fitter(numHiddenNodes=10)
costList = []

xRangeWide       = [0, 10]
numSamplesWide   = 50
batchSizeWide    = 50
trainDataWide    = DataSet(numSamplesWide, xRangeWide)
trainLoaderWide  = torch.utils.data.DataLoader(dataset=trainDataWide, batch_size=batchSizeWide, shuffle=True)

lossFn      = torch.nn.MSELoss()
optimiser   = torch.optim.Adam(network.parameters(), lr=1e-3)
epoch       = 0
numEpochs   = 1000
totalEpochs = 20000

start = time.time()
while epoch < totalEpochs:
    costList.extend(train(network, trainLoaderWide, lossFn, optimiser, numEpochs))
    epoch += numEpochs
end = time.time()


print("time = ", end-start, " s")
print(epoch, "epochs total, final cost = ", costList[-1])

plotNetwork(network, epoch)
plt.semilogy(costList)
plt.xlabel("Epochs",fontsize = 16)
plt.ylabel("Cost",fontsize = 16)
plt.title("Training Cost",fontsize = 16)
plt.show()

time =  25.50853967666626  s
20000 epochs total, final cost =  0.00015506161
MSE между аналитическим и полученным решениями =  1.049417733156588e-05