利用神经网络求解微分方程

最新推荐文章于 2024-05-15 15:26:03 发布

Rauchy

最新推荐文章于 2024-05-15 15:26:03 发布

阅读量1.9w

点赞数 9

分类专栏：数值优化与数学分析神经网络

本文链接：https://blog.csdn.net/wc13197389627/article/details/102644822

版权

神经网络同时被 2 个专栏收录

11 篇文章 3 订阅

订阅专栏

数值优化与数学分析

2 篇文章 2 订阅

订阅专栏

前言

最近在学习数值优化方面的问题，参考的资料为《数值最优化方法》，北大出版社的，高立编著，这是我们班的一个大佬推荐的。其中有提到用神经网络求解微分方程的为题，而且课后习题也有和神经网络算法相关的练习，于是打算提前研究一下，还可以锻炼自己的编程能力。

离散近似解法

最近在准备数学竞赛，刚好把微分方程又复习了一波。微分方程是由函数以及其导数组成的等式，一般而言，可以分为常微分方程(ODE)和偏微分方程(PDE)，常微分方程按照最高阶导数的阶数可以分为一阶，二阶甚至更高阶，按照函数及其导数的次数又可分为线性微分方程和非线性微分方程。求解微分方程一般有分离变量法和常数变易法等，但是这些方法可以求解的微分方程非常有限，对于更复杂的微分方程，求解解析解几乎是不可能的，于是在实际应用一般采用近似解法。

近似解法一般采用离散化的手段，化微分为差分，达到近似计算的效果。微分我们都很熟悉，函数f(x)在 $x_0$ 处的导数的定义如下，
$f(x_0)'=\lim_{x\to x_0} \frac{f(x)-f(x_0)}{x-x_0}=\lim _{\Delta x\to 0} \frac{\Delta f(x)}{\Delta x}$
当 $\Delta x$ 很小时，有，
$f(x)'\approx \frac{f(x+\Delta x)-f(x)}{\Delta x}$
上式也称为 $f (x)$ 一阶差商。利用类似的办法，我们可以得到高阶差商，偏导数也有类似的差商表示。高数里头我们学过泰勒公式，对于一个函数，我们可以利用它的导数将其展开，而且还可以估计泰勒展开式与原函数的误差。

考虑一个简单的一阶非线性微分方程 $f (x)^{'} + P (x) f (x) = Q (x)$
这类微分方程可以用常数变易法求解，这里我们采用近似的解法。简单起见，我们将x的取值范围限定为[0,1]，边界条件为f(0)=0(这里可以为其他值，为了便于叙述，我们假定为0)。求解过程如下：

我们将区间范围划分为n+1个格点，得到一个序列， $X=[x_0,x_1,x_2,...,x_n]$ ，其中 $x_0=0$ , $x_n=1$ ，分割精度为 $\Delta x=\frac{1}{n}$
利用中值定理，有 $f(x+\Delta x)\approx \Delta x f(x)'$ ，那么， $f(x_i)=\Delta x f(x_i)'$
这里我们利用前一点的导数值代替区间内的导数进行近似计算。
根据微分方程，有 $f (x)^{'} = Q (x) - P (x) f (x)$ 。根据初值条件 $f(x_0)=0$ , 进行一下递推计算
$f(x_{i-1})'=Q(x_{i-1})-P(x_{i-1})f(x_{i-1})$ $f(x_i)=f(x_{i-1})+\Delta xf(x_{i-1})'$

这样我们就完成了该微分方程在[0,1]上的近似计算。

实战一下，参考论文里给的一个例子是这样的，
$\frac{d}{dx}\Psi(x)+(x+\frac{1+3x^2}{1+x+x^3})\Psi(x)=x^3+2x+x^2\frac{1+3x^2}{1+x+x^3}$
$\Psi(0)=1$

这里参考的是medium里头作者的代码，如下：

import autograd.numpy as np
from autograd import grad 
import autograd.numpy.random as npr
from autograd.core import primitive
from matplotlib import pyplot as plt
%matplotlib inline

nx = 20
dx = 1. / nx

def A(x):
  '''
      Left part of initial equation
  '''
  return x + (1. + 3.*x**2) / (1. + x + x**3)


def B(x):
  '''
      Right part of initial equation
  '''
  return x**3 + 2.*x + x**2 * ((1. + 3.*x**2) / (1. + x + x**3))


def f(x, psy):
  '''
      d(psy)/dx = f(x, psy)
      This is f() function on the right
  '''
  return B(x) - psy * A(x)


def psy_analytic(x):
  '''
      Analytical solution of current problem
  '''
  return (np.exp((-x**2)/2.)) / (1. + x + x**3) + x**2

x_space = np.linspace(0, 1, nx)    
y_space = psy_analytic(x_space)
psy_fd = np.zeros_like(y_space)
psy_fd[0] = 1. # IC

for i in range(1, len(x_space)):
  psy_fd[i] = psy_fd[i-1] + B(x_space[i]) * dx - psy_fd[i-1] * A(x_space[i]) * dx
plt.figure()
analyticSolution, = plt.plot(x_space, y_space,'r-')
differentialSolution, = plt.plot(x_space, psy_fd,'g-')
plt.legend([analyticSolution,differentialSolution],['analytic solution','differential solution'],loc='upper left')
plt.show()

当格点的数量为10时，效果如下：
n=0
当格点的数量为20时，效果如下：
n=20
当区间分割得越细，计算的结果相较于解析解就越接近，当然，计算所需要的时间也随之线性增加。

对于更高阶的微分方程，就需要利用更高阶的差分了。

神经网络求解微分方程

接下来就是本篇文章的主题了，利用神经网络求解微分方程。神经网络从数学意义上来理解就是一个函数逼近器，它可以利用优化算法(back propagation)来逼近任意的函数，模型越复杂，表征能力越强，优化难度也越高。利用神经网络求解微分方程也是利用了神经网络模型的这一特点。

考虑这样一个微分方程，
$G(x,\Psi(x),\nabla\Psi(x),\nabla^2\Psi(x))=0,\forall x\in D$
其中， $x\in R^n,D\subset R^n$ 。边界条件域为 $\partial D$ 。为了满足边界条件，方程的可行解表示如下：
$\Psi_t(x,p)=\hat\Psi(x)+F(x)N(x,p)$
其中， $\Psi(x)=\hat\Psi(x), F(x)=0,\forall x\in \partial D$ ， $N (x, p)$ 是前馈神经网络，p是神经网络的参数。这个表达式很巧妙，因为它已经考虑了边界条件，所以我们在优化神经网络参数的时候就不用额外的考虑边界条件了。构造 $\Psi(x)$ 和 $F (x)$ 需要额外考虑，不过一般是采用关于x的多项式，后面我们会结合例子来进行讲解。

有了可行表达式，接下来就需要寻找优化的目标函数了。我们仍然将可行域 $D$ 进行离散化处理，将其划分为一个格点集 $\hat D=\{x^{(i)}\in D,i=1,2,...,m\}$ ，微分表达式在 $\hat D$ 上成立，即：
$G(x^{(i)},\Psi(x^{(i)}),\nabla\Psi(x^{(i)}),\nabla^2\Psi(x^{(i)}))=0,\forall x^{(i)}\in \hat D$
优化的目标函数为：
$J(p)=\sum_{i=1}^{m}G(x^{(i)},\Psi(x^{(i)}),\nabla\Psi(x^{(i)}),\nabla^2\Psi(x^{(i)}))^2$
让目标函数最小，也就是让可行解在格点集上无限接近微分表达式，最终，有
$p^*=\arg \min J(p)$

接下来就直接上实例了。优化目标函数需要进行微分运算，对于简单的函数，可以手动求微分，对于复杂的函数，手动求解积分较为复杂，一般采用Autograd来进行求解。实例中的微分表达式和刚才离散方法中的相同。优化的方法是采用梯度下降，同时控制下降速率。用到的是一个简单的三层神经网络：
神经网络结构
对于上述微分式，采用神经网络的计算方法得到的一个可行解表示如下：
$\Psi(x)=1+xN(x,p)$
这样，对 $\Psi(x)$ 求导如下：
$\frac{d}{dx}\Psi(x)=N(x,p)+x\frac{d}{dx}N(x,p)$
神经网络 $N (x, p)$ 的输入输出计算如下：

def neural_network(W, x):
    a1 = sigmoid(np.dot(x, W[0]))
    return np.dot(a1, W[1])

其中，x是输入，w是系数。 $N (x, p) = s i g m o i d (x W [0]) W [1]$

神经网络关于x的导函数 $\frac{d}{dx}N(x,p)$ 计算如下：

def d_neural_network_dx(W, x, k=1):
    return np.dot(np.dot(W[1].T, W[0].T**k), sigmoid_grad(x))

$\frac{d}{dx}N(x,p)=W[1]^TW[0]^T(sigmoid(x))$

损失函数计算如下:

def loss_function(W, x):
    loss_sum = 0.
    for xi in x:
        net_out = neural_network(W, xi)[0][0]
        psy_t = 1. + xi * net_out
        d_net_out = d_neural_network_dx(W, xi)[0][0]
        d_psy_t = net_out + xi * d_net_out
        func = f(xi, psy_t)       
        err_sqr = (d_psy_t - func)**2
        loss_sum += err_sqr
    return loss_sum

这里的计算稍复杂一些，需要指出的是，上式中的 $d\_psy\_t-func$ 就是我们之前提到的微分表达式G，有兴趣的话可以查看完整的代码，为了缩短文章的篇幅，我没有把全部的代码贴上来。
最后运行的效果如下：
n=100
这是n=100时运行的效果，由于区间为[0,1]，此时利用离散的解法可以很好的近似解，而利用神经网络的解法不但时间长，而且效果没有近似解法好。