大家好，欢迎来到IT知识分享网。
非线性薛定谔方程的形式
\(ih_t + 0.5h_{xx} + |h|^2 h = 0, \ \ x \in [-5,5],t\in [0,\pi/2]\)
\(h(0,x) = 2sech(x)\)
\(h(t,-5) = 2sech(x)\)
\(h_x(t,-5) = h_x(t,5)\)
1. 算法思路
\(MSE = MSE_0 + MSE_b + MSE_f\)
其中，
\(MSE_0 = \frac{1}{N_0} \sum\limits_{i=1}^{N_0}|h(0,x_0^i)-h_0^i|^2···(1)\)
\(MSE_b = \frac{1}{N_b}\sum\limits_{i=1}^{N_b}|h^i(t_b^i,-5)-h^i(t_b^i,5)|^2 + |h_x^i(t_b^i,-5)-h_x^i(t_b^i,5)|^2···(2)\)
\(MSE_f = \frac{1}{N_f}\sum\limits_{i=1}^{N_f}|f(t_f^i,x_f^i)-0|^2···(3)\)
损失函数分为三部分：
- 输入（x，0），MSE_0表示预测出的初始值与真实初值之间的差别
- 输入(x边界，随机t)，MSE_b表示预测值满足周期性边界条件的程度
- 在定义域内随机选取时空点（x,t），代入神经网络得到的预测值与将（x,t）直接代入方程之间的误差记为MSE_f
2. 代码实现
2.1 简要解读
主干部分
以下是建立神经网络的主要代码：
# 建立神经网络图
self.u0_pred, self.v0_pred, _ , _ = self.net_uv(self.x0_tf, self.t0_tf)
self.u_lb_pred, self.v_lb_pred, self.u_x_lb_pred, self.v_x_lb_pred = self.net_uv(self.x_lb_tf, self.t_lb_tf)
self.u_ub_pred, self.v_ub_pred, self.u_x_ub_pred, self.v_x_ub_pred = self.net_uv(self.x_ub_tf, self.t_ub_tf)
self.f_u_pred, self.f_v_pred = self.net_f_uv(self.x_f_tf, self.t_f_tf)
建立神经网络个神经元之间的计算关系：
def neural_net(self, X, weights, biases):
num_layers = len(weights) + 1
H = 2.0*(X - self.lb)/(self.ub - self.lb) - 1.0
for l in range(0,num_layers-2):
W = weights[l]
b = biases[l]
H = tf.tanh(tf.add(tf.matmul(H, W), b))
W = weights[-1]
b = biases[-1]
Y = tf.add(tf.matmul(H, W), b)
return Y
计算初边值点：
def net_uv(self, x, t):
X = tf.concat([x,t],1)
uv = self.neural_net(X, self.weights, self.biases)
u = uv[:,0:1]
v = uv[:,1:2]
u_x = tf.gradients(u, x)[0]
v_x = tf.gradients(v, x)[0]
return u, v, u_x, v_x
计算内部配置点：
def net_f_uv(self, x, t):
u, v, u_x, v_x = self.net_uv(x,t)
u_t = tf.gradients(u, t)[0]
u_xx = tf.gradients(u_x, x)[0]
v_t = tf.gradients(v, t)[0]
v_xx = tf.gradients(v_x, x)[0]
f_u = u_t + 0.5*v_xx + (u**2 + v**2)*v
f_v = v_t - 0.5*u_xx - (u**2 + v**2)*u
return f_u, f_v
LOSS函数按照公式定义如下：
# Loss
self.loss = tf.reduce_mean(tf.square(self.u0_tf - self.u0_pred)) + \
tf.reduce_mean(tf.square(self.v0_tf - self.v0_pred)) + \
tf.reduce_mean(tf.square(self.u_lb_pred - self.u_ub_pred)) + \
tf.reduce_mean(tf.square(self.v_lb_pred - self.v_ub_pred)) + \
tf.reduce_mean(tf.square(self.u_x_lb_pred - self.u_x_ub_pred)) + \
tf.reduce_mean(tf.square(self.v_x_lb_pred - self.v_x_ub_pred)) + \
tf.reduce_mean(tf.square(self.f_u_pred)) + \
tf.reduce_mean(tf.square(self.f_v_pred))
使用了两个优化器，定义如下：
# Optimizers
self.optimizer = tf.contrib.opt.ScipyOptimizerInterface(self.loss,
method = 'L-BFGS-B',
options = {'maxiter': 50000,
'maxfun': 50000,
'maxcor': 50,
'maxls': 50,
'ftol' : 1.0 * np.finfo(float).eps})
self.optimizer_Adam = tf.train.AdamOptimizer()
self.train_op_Adam = self.optimizer_Adam.minimize(self.loss)
在函数train中使用如下：
def train(self, nIter): #nIter是训练次数，会在main函数中给出
# 字典
tf_dict = {self.x0_tf: self.x0, self.t0_tf: self.t0,
self.u0_tf: self.u0, self.v0_tf: self.v0,
self.x_lb_tf: self.x_lb, self.t_lb_tf: self.t_lb,
self.x_ub_tf: self.x_ub, self.t_ub_tf: self.t_ub,
self.x_f_tf: self.x_f, self.t_f_tf: self.t_f}
start_time = time.time()
for it in range(nIter):
self.sess.run(self.train_op_Adam, tf_dict) #使用Adam优化器预处理，因为其收敛快，但是精度较低
# Print
if it % 10 == 0:
elapsed = time.time() - start_time
loss_value = self.sess.run(self.loss, tf_dict)
print('It: %d, Loss: %.3e, Time: %.2f' %
(it, loss_value, elapsed))
start_time = time.time()
# 使用L-BFGS-B优化器进一步使Loss减小，因为其收敛较慢，但是精度较高
self.optimizer.minimize(self.sess,
feed_dict = tf_dict,
fetches = [self.loss],
loss_callback = self.callback)
初始化
权重和偏置矩阵的初始化如下：
def initialize_NN(self, layers):
weights = []
biases = []
num_layers = len(layers)
for l in range(0,num_layers-1):
W = self.xavier_init(size=[layers[l], layers[l+1]]) # 权重使用Xavier方法初始化
b = tf.Variable(tf.zeros([1,layers[l+1]], dtype=tf.float32), dtype=tf.float32) # 偏置直接初始化为零
weights.append(W)
biases.append(b)
return weights, biases
def xavier_init(self, size):
in_dim = size[0]
out_dim = size[1]
xavier_stddev = np.sqrt(2/(in_dim + out_dim))
return tf.Variable(tf.truncated_normal([in_dim, out_dim], stddev=xavier_stddev), dtype=tf.float32)
测试
main函数中的测试代码部分：
# 用于构建神经网络的只是部分数据，而这里的X_star是所有的（x,t）
u_pred, v_pred, f_u_pred, f_v_pred = model.predict(X_star) # 代入神经网络得到预测值
h_pred = np.sqrt(u_pred**2 + v_pred**2) # 复数 = sqrt(实部^2 + 虚部^2)
# 衡量预测值与真实值之间的误差
error_u = np.linalg.norm(u_star-u_pred,2)/np.linalg.norm(u_star,2)
error_v = np.linalg.norm(v_star-v_pred,2)/np.linalg.norm(v_star,2)
error_h = np.linalg.norm(h_star-h_pred,2)/np.linalg.norm(h_star,2)
# 输出误差
print('Error u: %e' % (error_u))
print('Error v: %e' % (error_v))
print('Error h: %e' % (error_h))
class PhysicsInformedNN中的相关部分：
def predict(self, X_star):
tf_dict = {self.x0_tf: X_star[:,0:1], self.t0_tf: X_star[:,1:2]}
u_star = self.sess.run(self.u0_pred, tf_dict)
v_star = self.sess.run(self.v0_pred, tf_dict)
tf_dict = {self.x_f_tf: X_star[:,0:1], self.t_f_tf: X_star[:,1:2]}
f_u_star = self.sess.run(self.f_u_pred, tf_dict)
f_v_star = self.sess.run(self.f_v_pred, tf_dict)
return u_star, v_star, f_u_star, f_v_star
完整代码按顺序解读
点击查看代码
"""
@author: Maziar Raissi
"""
import sys
sys.path.insert(0, '../../Utilities/')
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
from scipy.interpolate import griddata
from pyDOE import lhs
from plotting import newfig, savefig
from mpl_toolkits.mplot3d import Axes3D
import time
import matplotlib.gridspec as gridspec
from mpl_toolkits.axes_grid1 import make_axes_locatable
np.random.seed(1234)
tf.set_random_seed(1234)
class PhysicsInformedNN:
# Initialize the class
def __init__(self, x0, u0, v0, tb, X_f, layers, lb, ub):
X0 = np.concatenate((x0, 0*x0), 1) # (x0, 0)
X_lb = np.concatenate((0*tb + lb[0], tb), 1) # (lb[0], tb)
X_ub = np.concatenate((0*tb + ub[0], tb), 1) # (ub[0], tb)
self.lb = lb
self.ub = ub
self.x0 = X0[:,0:1]
self.t0 = X0[:,1:2]
self.x_lb = X_lb[:,0:1]
self.t_lb = X_lb[:,1:2]
self.x_ub = X_ub[:,0:1]
self.t_ub = X_ub[:,1:2]
self.x_f = X_f[:,0:1]
self.t_f = X_f[:,1:2]
self.u0 = u0
self.v0 = v0
# Initialize NNs
self.layers = layers
self.weights, self.biases = self.initialize_NN(layers)
# tf Placeholders
self.x0_tf = tf.placeholder(tf.float32, shape=[None, self.x0.shape[1]])
self.t0_tf = tf.placeholder(tf.float32, shape=[None, self.t0.shape[1]])
self.u0_tf = tf.placeholder(tf.float32, shape=[None, self.u0.shape[1]])
self.v0_tf = tf.placeholder(tf.float32, shape=[None, self.v0.shape[1]])
self.x_lb_tf = tf.placeholder(tf.float32, shape=[None, self.x_lb.shape[1]])
self.t_lb_tf = tf.placeholder(tf.float32, shape=[None, self.t_lb.shape[1]])
self.x_ub_tf = tf.placeholder(tf.float32, shape=[None, self.x_ub.shape[1]])
self.t_ub_tf = tf.placeholder(tf.float32, shape=[None, self.t_ub.shape[1]])
self.x_f_tf = tf.placeholder(tf.float32, shape=[None, self.x_f.shape[1]])
self.t_f_tf = tf.placeholder(tf.float32, shape=[None, self.t_f.shape[1]])
# tf Graphs
self.u0_pred, self.v0_pred, _ , _ = self.net_uv(self.x0_tf, self.t0_tf)
self.u_lb_pred, self.v_lb_pred, self.u_x_lb_pred, self.v_x_lb_pred = self.net_uv(self.x_lb_tf, self.t_lb_tf)
self.u_ub_pred, self.v_ub_pred, self.u_x_ub_pred, self.v_x_ub_pred = self.net_uv(self.x_ub_tf, self.t_ub_tf)
self.f_u_pred, self.f_v_pred = self.net_f_uv(self.x_f_tf, self.t_f_tf)
# Loss
self.loss = tf.reduce_mean(tf.square(self.u0_tf - self.u0_pred)) + \
tf.reduce_mean(tf.square(self.v0_tf - self.v0_pred)) + \
tf.reduce_mean(tf.square(self.u_lb_pred - self.u_ub_pred)) + \
tf.reduce_mean(tf.square(self.v_lb_pred - self.v_ub_pred)) + \
tf.reduce_mean(tf.square(self.u_x_lb_pred - self.u_x_ub_pred)) + \
tf.reduce_mean(tf.square(self.v_x_lb_pred - self.v_x_ub_pred)) + \
tf.reduce_mean(tf.square(self.f_u_pred)) + \
tf.reduce_mean(tf.square(self.f_v_pred))
# Optimizers
self.optimizer = tf.contrib.opt.ScipyOptimizerInterface(self.loss,
method = 'L-BFGS-B',
options = {'maxiter': 50000,
'maxfun': 50000,
'maxcor': 50,
'maxls': 50,
'ftol' : 1.0 * np.finfo(float).eps})
self.optimizer_Adam = tf.train.AdamOptimizer()
self.train_op_Adam = self.optimizer_Adam.minimize(self.loss)
# tf session
self.sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True,
log_device_placement=True))
init = tf.global_variables_initializer()
self.sess.run(init)
def initialize_NN(self, layers):
weights = []
biases = []
num_layers = len(layers)
for l in range(0,num_layers-1):
W = self.xavier_init(size=[layers[l], layers[l+1]])
b = tf.Variable(tf.zeros([1,layers[l+1]], dtype=tf.float32), dtype=tf.float32)
weights.append(W)
biases.append(b)
return weights, biases
def xavier_init(self, size):
in_dim = size[0]
out_dim = size[1]
xavier_stddev = np.sqrt(2/(in_dim + out_dim))
return tf.Variable(tf.truncated_normal([in_dim, out_dim], stddev=xavier_stddev), dtype=tf.float32)
def neural_net(self, X, weights, biases):
num_layers = len(weights) + 1
H = 2.0*(X - self.lb)/(self.ub - self.lb) - 1.0
for l in range(0,num_layers-2):
W = weights[l]
b = biases[l]
H = tf.tanh(tf.add(tf.matmul(H, W), b))
W = weights[-1]
b = biases[-1]
Y = tf.add(tf.matmul(H, W), b)
return Y
def net_uv(self, x, t):
X = tf.concat([x,t],1)
uv = self.neural_net(X, self.weights, self.biases)
u = uv[:,0:1]
v = uv[:,1:2]
u_x = tf.gradients(u, x)[0]
v_x = tf.gradients(v, x)[0]
return u, v, u_x, v_x
def net_f_uv(self, x, t):
u, v, u_x, v_x = self.net_uv(x,t)
u_t = tf.gradients(u, t)[0]
u_xx = tf.gradients(u_x, x)[0]
v_t = tf.gradients(v, t)[0]
v_xx = tf.gradients(v_x, x)[0]
# 代入公式
f_u = u_t + 0.5*v_xx + (u**2 + v**2)*v # 因为有虚数i,因此f_u的第三项：原本的v乘以i之后就变成了实部的一部分，因此这里需要是(u**2 + v**2)*v
f_v = v_t - 0.5*u_xx - (u**2 + v**2)*u
return f_u, f_v
def callback(self, loss):
print('Loss:', loss)
def train(self, nIter):
tf_dict = {self.x0_tf: self.x0, self.t0_tf: self.t0,
self.u0_tf: self.u0, self.v0_tf: self.v0,
self.x_lb_tf: self.x_lb, self.t_lb_tf: self.t_lb,
self.x_ub_tf: self.x_ub, self.t_ub_tf: self.t_ub,
self.x_f_tf: self.x_f, self.t_f_tf: self.t_f}
start_time = time.time()
for it in range(nIter):
self.sess.run(self.train_op_Adam, tf_dict)
# Print
if it % 10 == 0:
elapsed = time.time() - start_time
loss_value = self.sess.run(self.loss, tf_dict)
print('It: %d, Loss: %.3e, Time: %.2f' %
(it, loss_value, elapsed))
start_time = time.time()
self.optimizer.minimize(self.sess,
feed_dict = tf_dict,
fetches = [self.loss],
loss_callback = self.callback)
# 用于测试神经网络性能
def predict(self, X_star):
tf_dict = {self.x0_tf: X_star[:,0:1], self.t0_tf: X_star[:,1:2]}
u_star = self.sess.run(self.u0_pred, tf_dict)
v_star = self.sess.run(self.v0_pred, tf_dict)
tf_dict = {self.x_f_tf: X_star[:,0:1], self.t_f_tf: X_star[:,1:2]}
f_u_star = self.sess.run(self.f_u_pred, tf_dict)
f_v_star = self.sess.run(self.f_v_pred, tf_dict)
return u_star, v_star, f_u_star, f_v_star
if __name__ == "__main__": # python知识，加上这句话，单独运行这个文件的时候会运行main，也可以直接调用函数PhysicsInformedNN，而无需管main函数
noise = 0.0
# Doman bounds
lb = np.array([-5.0, 0.0])
ub = np.array([5.0, np.pi/2])
N0 = 50 #初始点
N_b = 50 #边界点
N_f = 20000 #内部配置点
layers = [2, 100, 100, 100, 100, 2]
data = scipy.io.loadmat('../Data/NLS.mat')
t = data['tt'].flatten()[:,None]
x = data['x'].flatten()[:,None]
Exact = data['uu']
Exact_u = np.real(Exact)
Exact_v = np.imag(Exact)
Exact_h = np.sqrt(Exact_u**2 + Exact_v**2)
X, T = np.meshgrid(x,t)
X_star = np.hstack((X.flatten()[:,None], T.flatten()[:,None]))
u_star = Exact_u.T.flatten()[:,None]
v_star = Exact_v.T.flatten()[:,None]
h_star = Exact_h.T.flatten()[:,None]
###########################
# 从数据x中随机抽取N0=50个点
idx_x = np.random.choice(x.shape[0], N0, replace=False) # 指标index
x0 = x[idx_x,:] # 从数据x中随机抽取50个点，命名为x0
u0 = Exact_u[idx_x,0:1] # 实部
v0 = Exact_v[idx_x,0:1] # 虚部
# 从数据t中随机抽取N_b=50个点命名为tb(与x0并不一一对应)
idx_t = np.random.choice(t.shape[0], N_b, replace=False)
tb = t[idx_t,:]
# 从上下界中随机抽取N_f=20000个点，命名为X_f
X_f = lb + (ub-lb)*lhs(2, N_f) # lhs:拉丁超立方抽样
model = PhysicsInformedNN(x0, u0, v0, tb, X_f, layers, lb, ub)
start_time = time.time()
model.train(50000)
elapsed = time.time() - start_time
print('Training time: %.4f' % (elapsed))
u_pred, v_pred, f_u_pred, f_v_pred = model.predict(X_star)
h_pred = np.sqrt(u_pred**2 + v_pred**2)
error_u = np.linalg.norm(u_star-u_pred,2)/np.linalg.norm(u_star,2)
error_v = np.linalg.norm(v_star-v_pred,2)/np.linalg.norm(v_star,2)
error_h = np.linalg.norm(h_star-h_pred,2)/np.linalg.norm(h_star,2)
print('Error u: %e' % (error_u))
print('Error v: %e' % (error_v))
print('Error h: %e' % (error_h))
U_pred = griddata(X_star, u_pred.flatten(), (X, T), method='cubic')
V_pred = griddata(X_star, v_pred.flatten(), (X, T), method='cubic')
H_pred = griddata(X_star, h_pred.flatten(), (X, T), method='cubic')
FU_pred = griddata(X_star, f_u_pred.flatten(), (X, T), method='cubic')
FV_pred = griddata(X_star, f_v_pred.flatten(), (X, T), method='cubic')
######################################################################
############################# Plotting ###############################
######################################################################
# matplotlib还没学会，待更......
X0 = np.concatenate((x0, 0*x0), 1) # (x0, 0)
X_lb = np.concatenate((0*tb + lb[0], tb), 1) # (lb[0], tb)
X_ub = np.concatenate((0*tb + ub[0], tb), 1) # (ub[0], tb)
X_u_train = np.vstack([X0, X_lb, X_ub])
fig, ax = newfig(1.0, 0.9)
ax.axis('off')
####### Row 0: h(t,x) ##################
gs0 = gridspec.GridSpec(1, 2)
gs0.update(top=1-0.06, bottom=1-1/3, left=0.15, right=0.85, wspace=0)
ax = plt.subplot(gs0[:, :])
h = ax.imshow(H_pred.T, interpolation='nearest', cmap='YlGnBu',
extent=[lb[1], ub[1], lb[0], ub[0]],
origin='lower', aspect='auto')
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(h, cax=cax)
ax.plot(X_u_train[:,1], X_u_train[:,0], 'kx', label = 'Data (%d points)' % (X_u_train.shape[0]), markersize = 4, clip_on = False)
line = np.linspace(x.min(), x.max(), 2)[:,None]
ax.plot(t[75]*np.ones((2,1)), line, 'k--', linewidth = 1)
ax.plot(t[100]*np.ones((2,1)), line, 'k--', linewidth = 1)
ax.plot(t[125]*np.ones((2,1)), line, 'k--', linewidth = 1)
ax.set_xlabel('$t$')
ax.set_ylabel('$x$')
leg = ax.legend(frameon=False, loc = 'best')
# plt.setp(leg.get_texts(), color='w')
ax.set_title('$|h(t,x)|$', fontsize = 10)
####### Row 1: h(t,x) slices ##################
gs1 = gridspec.GridSpec(1, 3)
gs1.update(top=1-1/3, bottom=0, left=0.1, right=0.9, wspace=0.5)
ax = plt.subplot(gs1[0, 0])
ax.plot(x,Exact_h[:,75], 'b-', linewidth = 2, label = 'Exact')
ax.plot(x,H_pred[75,:], 'r--', linewidth = 2, label = 'Prediction')
ax.set_xlabel('$x$')
ax.set_ylabel('$|h(t,x)|$')
ax.set_title('$t = %.2f$' % (t[75]), fontsize = 10)
ax.axis('square')
ax.set_xlim([-5.1,5.1])
ax.set_ylim([-0.1,5.1])
ax = plt.subplot(gs1[0, 1])
ax.plot(x,Exact_h[:,100], 'b-', linewidth = 2, label = 'Exact')
ax.plot(x,H_pred[100,:], 'r--', linewidth = 2, label = 'Prediction')
ax.set_xlabel('$x$')
ax.set_ylabel('$|h(t,x)|$')
ax.axis('square')
ax.set_xlim([-5.1,5.1])
ax.set_ylim([-0.1,5.1])
ax.set_title('$t = %.2f$' % (t[100]), fontsize = 10)
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.8), ncol=5, frameon=False)
ax = plt.subplot(gs1[0, 2])
ax.plot(x,Exact_h[:,125], 'b-', linewidth = 2, label = 'Exact')
ax.plot(x,H_pred[125,:], 'r--', linewidth = 2, label = 'Prediction')
ax.set_xlabel('$x$')
ax.set_ylabel('$|h(t,x)|$')
ax.axis('square')
ax.set_xlim([-5.1,5.1])
ax.set_ylim([-0.1,5.1])
ax.set_title('$t = %.2f$' % (t[125]), fontsize = 10)
# savefig('./figures/NLS')
免责声明：本站所有文章内容,图片，视频等均是来源于用户投稿和互联网及文摘转载整编而成，不代表本站观点，不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益，请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权，请联系删除，如若转载，请注明出处：https://yundeesoft.com/30508.html
