当前位置:   article > 正文

Python 机器学习求解 PDE 学习项目——PINN 求解二维 Poisson 方程

Python 机器学习求解 PDE 学习项目——PINN 求解二维 Poisson 方程

本文使用 TensorFlow 1.15 环境搭建深度神经网络(PINN)求解二维 Poisson 方程:

模型问题

− Δ u = f in  Ω , u = g on  Γ : = ∂ Ω . Δu=fin Ω,u=gon Γ:=Ω. Δuu=f=gin Ω,on Γ:=Ω.
其中 Ω = [ X a , X b ] × [ Y a , Y b ] \Omega = [X_a,X_b]\times[Y_a,Y_b] Ω=[Xa,Xb]×[Ya,Yb] 是一个二维矩形区域, Δ u = u x x + u y y , g \Delta u = u_{xx}+u_{yy}, g Δu=uxx+uyy,g 是边界条件给定的函数,可以非零.

在这里插入图片描述

代码展现

二维PINN 与一维的整体框架是类似的,只是数据的维度升高了,为了读者方便这里直接展示完整代码,每段代码都添加了注释帮助理解:

检查 TF 版本号:

# PINN 求解 2D Poisson 方程
import tensorflow as tf
print(tf.__version__)
  • 1
  • 2
  • 3

主要函数

import os
#tensorflow-intel automatically set the TF_ENABLE_ONEDNN_OPTS=1
#os.environ['TF_ENABLE_ONEDNN_OPTS'] = '0'
# Here, TF_ENABLE_ONEDNN_OPTS=0 should be above import tensorflow as tf
import tensorflow as tf
import numpy as np
import time
import matplotlib.pyplot as plt
import scipy.io
import math
# 定义数据集类,用于生成训练所需的数据
class Dataset:
    def __init__(self, x_range, y_range, N_res, N_bx, N_by, Nx, Ny, xa, xb, ya, yb):
        self.x_range = x_range  # x 轴范围
        self.y_range = y_range  # y 轴范围
        self.N_res = N_res  # 方程残差点数量
        self.N_bx = N_bx  # x 方向边界条件点数量
        self.N_by = N_by  # y 方向边界条件点数量
        self.Nx = Nx  # x 方向网格数量
        self.Ny = Ny  # y 方向网格数量
        self.xa = xa  # x 方向左边界
        self.xb = xb  # x 方向右边界
        self.ya = ya  # y 方向下边界
        self.yb = yb  # y 方向上边界

    # 定义边界条件函数
    # 可以求解非齐次 Dirichlet 边界条件
    def bc(self, X_b):
        U_bc = Exact(self.xa, self.xb, self.ya, self.yb)
        u_bc = U_bc.u_exact(X_b)
        return u_bc

    # 生成数据:残差点和边界条件点
    def build_data(self):
        x0, x1 = self.x_range
        y0, y1 = self.y_range

        Xmin = np.hstack((x0, y0))
        Xmax = np.hstack((x1, y1))

        ## 如果使用均匀网格,代码如下:
        """
        x_ = np.linspace(x0, x1, self.Nx).reshape((-1, 1))
        y_ = np.linspace(y0, y1, self.Ny).reshape((-1, 1))
        x, y = np.meshgrid(x_, y_)
        x = np.reshape(x, (-1, 1))
        y = np.reshape(y, (-1, 1))
        xy = np.hstack((x, y))
        X_res_input = xy
        """

        ## 为方程生成随机残差点
        x_res = x0 + (x1 - x0) * np.random.rand(self.N_res, 1)
        y_res = y0 + (y1 - y0) * np.random.rand(self.N_res, 1)
        X_res_input = np.hstack((x_res, y_res))

        # 生成 x = xa, xb 的边界条件点
        y_b = y0 + (y1 - y0) * np.random.rand(self.N_by, 1)
        x_b0 = x0 * np.ones_like(y_b)
        x_b1 = x1 * np.ones_like(y_b)
        X_b0_input = np.hstack((x_b0, y_b))
        X_b1_input = np.hstack((x_b1, y_b))

        # 生成 y = ya, yb 的边界条件点
        x_b = x0 + (x1 - x0) * np.random.rand(self.N_bx, 1)
        y_b0 = y0 * np.ones_like(x_b)
        y_b1 = y1 * np.ones_like(x_b)
        Y_b0_input = np.hstack((x_b, y_b0))
        Y_b1_input = np.hstack((x_b, y_b1))
        
        return X_res_input, X_b0_input, X_b1_input, Y_b0_input, Y_b1_input, Xmin, Xmax

# 定义精确解类,用于计算精确解
class Exact:
    def __init__(self, xa, xb, ya, yb):
        self.xa = xa  # x 方向左边界
        self.xb = xb  # x 方向右边界
        self.ya = ya  # y 方向下边界
        self.yb = yb  # y 方向上边界

    # 精确解函数
    def u_exact(self, X):
        x = X[:, 0:1]
        y = X[:, 1:2]

        u = np.sin(2 * np.pi * x  ) * np.sin(2 * np.pi * y )
        return u

class Train:
    def __init__(self, train_dict):
        self.train_dict = train_dict  # 训练数据
        self.step = 0  # 训练步数

    # 打印训练损失
    def callback(self, loss_value):
        self.step += 1
        if self.step % 200 == 0:
            print(f'Loss: {loss_value:.4e}')

    # 使用 Adam 和 L-BFGS 优化器进行训练
    def nntrain(self, sess, u_pred, loss, train_adam, train_lbfgs):
        n = 0
        max_steps = 1000
        loss_threshold = 4.0e-4
        current_loss = 1.0
        while n < max_steps and current_loss > loss_threshold:
            n += 1
            u_, current_loss, _ = sess.run([u_pred, loss, train_adam], feed_dict=self.train_dict)
            # 每2^n步打印一次损失并绘制结果  
            if math.isclose(math.fmod(math.log2(n), 1), 0, abs_tol=1e-9): 
                print(f'Steps: {n}, loss: {current_loss:.4e}')

        train_lbfgs.minimize(sess, feed_dict=self.train_dict, fetches=[loss], loss_callback=self.callback)

class DNN:
    def __init__(self, layer_sizes, Xmin, Xmax):
        self.layer_sizes = layer_sizes  # 每层的节点数
        self.Xmin = Xmin  # 输入范围最小值
        self.Xmax = Xmax  # 输入范围最大值
    
    # 初始化神经网络的权重和偏置
    def hyper_initial(self):
        num_layers = len(self.layer_sizes)
        weights = []
        biases = []
        for l in range(1, num_layers):
            in_dim = self.layer_sizes[l-1]
            out_dim = self.layer_sizes[l]
            std = np.sqrt(2 / (in_dim + out_dim))
            weight = tf.Variable(tf.random_normal(shape=[in_dim, out_dim], stddev=std))
            bias = tf.Variable(tf.zeros(shape=[1, out_dim]))
            weights.append(weight)
            biases.append(bias)
        return weights, biases

    # 构建前馈神经网络
    def fnn(self, X, weights, biases):
        A = 2.0 * (X - self.Xmin) / (self.Xmax - self.Xmin) - 1.0  # 归一化输入
        num_layers = len(weights)
        for i in range(num_layers - 1):
            A = tf.tanh(tf.add(tf.matmul(A, weights[i]), biases[i]))  # 隐藏层激活函数
        Y = tf.add(tf.matmul(A, weights[-1]), biases[-1])  # 输出层
        return Y
    
    # 构建用于求解 Poisson 方程的神经网络
    def pdenn(self, x, y, weights, biases):
        u = self.fnn(tf.concat([x, y], 1), weights, biases)  # 前馈网络输出
        u_x = tf.gradients(u, x)[0]  # u 对 x 的一阶导数
        u_xx = tf.gradients(u_x, x)[0]  # u 对 x 的二阶导数
        u_y = tf.gradients(u, y)[0]  # u 对 y 的一阶导数
        u_yy = tf.gradients(u_y, y)[0]  # u 对 y 的二阶导数
        
        # 源项函数
        rhs_func = 8 * np.pi**2 * tf.sin(2 * np.pi * x  ) * tf.sin(2 * np.pi * y )

        # 残差项
        residual = -(u_xx + u_yy) - rhs_func

        return residual

def compute_errors(u_pred, u_exact):
    """
    计算数值解与精确解之间的 L2 误差和最大模误差
    :param u_pred: 数值解
    :param u_exact: 精确解
    :return: L2 误差和最大模误差
    """
    # 计算 L2 误差
    L2_error = np.sqrt(np.mean((u_pred - u_exact) ** 2))
    
    # 计算最大模误差
    max_error = np.max(np.abs(u_pred - u_exact))
    
    return L2_error, max_error

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129
  • 130
  • 131
  • 132
  • 133
  • 134
  • 135
  • 136
  • 137
  • 138
  • 139
  • 140
  • 141
  • 142
  • 143
  • 144
  • 145
  • 146
  • 147
  • 148
  • 149
  • 150
  • 151
  • 152
  • 153
  • 154
  • 155
  • 156
  • 157
  • 158
  • 159
  • 160
  • 161
  • 162
  • 163
  • 164
  • 165
  • 166
  • 167
  • 168
  • 169
  • 170
  • 171
  • 172
  • 173
  • 174
  • 175

画图以及保存图片:

# 检查保存路径是否存在,如果不存在则创建
save_path = './Output'
if not os.path.exists(save_path):
    os.makedirs(save_path)

# 定义保存和绘图类
class SavePlot:
    def __init__(self, session, x_range, y_range, num_x_points, num_y_points, xa, xb, ya, yb):
        self.x_range = x_range  # x 轴范围
        self.y_range = y_range  # y 轴范围
        self.num_x_points = num_x_points  # x 方向上的测试点数量
        self.num_y_points = num_y_points  # y 方向上的测试点数量
        self.session = session  # TensorFlow 会话
        self.xa = xa  # x 方向左边界
        self.xb = xb  # x 方向右边界
        self.ya = ya  # y 方向下边界
        self.yb = yb  # y 方向上边界

    # 保存并绘制预测和精确解
    def save_and_plot(self, u_pred, x_res_train, y_res_train):
        # 生成测试点
        x_test = np.linspace(self.x_range[0], self.x_range[1], self.num_x_points).reshape((-1, 1))
        y_test = np.linspace(self.y_range[0], self.y_range[1], self.num_y_points).reshape((-1, 1))
        x_test_grid, y_test_grid = np.meshgrid(x_test, y_test)
        x_test_grid = np.reshape(x_test_grid, (-1, 1))
        y_test_grid = np.reshape(y_test_grid, (-1, 1))
        
        # 创建测试字典
        test_feed_dict = {x_res_train: x_test_grid, y_res_train: y_test_grid}
        
        # 在测试网格上进行预测
        u_test = self.session.run(u_pred, feed_dict=test_feed_dict)
        u_test = np.reshape(u_test, (y_test.shape[0], x_test.shape[0]))
        u_test = np.transpose(u_test)

        # 保存预测结果到文件
        np.savetxt(os.path.join(save_path, 'u_pred.txt'), u_test, fmt='%e')

        # 绘制预测结果并保存图片
        plt.imshow(u_test, cmap='rainbow', aspect='auto')
        plt.colorbar()
        plt.title('Numerical Solution')
        plt.xlabel('X-axis')
        plt.ylabel('Y-axis')
        plt.savefig(os.path.join(save_path, 'u_pred.png'))
        plt.show()
        plt.close()

        # 计算并保存精确解
        exact_solution = Exact(self.xa, self.xb, self.ya, self.yb)
        u_exact = exact_solution.u_exact(np.hstack((x_test_grid, y_test_grid)))
        u_exact = np.reshape(u_exact, (y_test.shape[0], x_test.shape[0]))
        u_exact = np.transpose(u_exact)
        np.savetxt(os.path.join(save_path, 'u_exact.txt'), u_exact, fmt='%e')

        # 绘制精确解并保存图片
        plt.imshow(u_exact, cmap='rainbow', aspect='auto')
        plt.colorbar()
        plt.title('Exact Solution')
        plt.xlabel('X-axis')
        plt.ylabel('Y-axis')
        plt.savefig(os.path.join(save_path, 'u_exact.png'))
        plt.show()
        plt.close()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64

下面是主程序:

import os
import tensorflow as tf
import numpy as np
import time
import matplotlib.pyplot as plt

# 设置随机种子以确保可重复性
np.random.seed(1234)
tf.set_random_seed(1234)

def main():
    # 定义计算域范围
    x_range = [-0.5, 1.5]
    y_range = [-1.0, 1.0]
    
    # 网格点数量
    num_x_points = 101
    num_y_points = 101
    
    # 残差点和边界点数量
    num_residual_points = 8000
    num_boundary_x_points = 100
    num_boundary_y_points = 100
    
    # 边界范围
    xa = x_range[0]
    xb = x_range[1]
    ya = y_range[0]
    yb = y_range[1]

    # 创建数据集对象
    data = Dataset(x_range, y_range, num_residual_points, num_boundary_x_points, num_boundary_y_points, num_x_points, num_y_points, xa, xb, ya, yb)
    
    # 生成数据
    X_res, X_b0, X_b1, Y_b0, Y_b1, Xmin, Xmax = data.build_data()

    # 定义神经网络的层结构
    layers = [2] + 5 * [40] + [1]

    # 定义输入占位符
    x_res_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    y_res_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    X_x_b0_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    X_y_b0_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    X_x_b1_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    X_y_b1_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    Y_x_b0_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    Y_y_b0_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    Y_x_b1_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)
    Y_y_b1_train = tf.placeholder(shape=[None, 1], dtype=tf.float32)

    # 创建物理信息神经网络(PINN)
    pinn = DNN(layers, Xmin, Xmax)
    weights, biases = pinn.hyper_initial()

    # 预测解
    u_pred = pinn.fnn(tf.concat([x_res_train, y_res_train], 1), weights, biases)

    # 计算残差
    f_pred = pinn.pdenn(x_res_train, y_res_train, weights, biases)

    # 边界条件预测 (x = xa, xb)
    u_x_b0_pred = pinn.fnn(tf.concat([X_x_b0_train, X_y_b0_train], 1), weights, biases)
    u_x_b1_pred = pinn.fnn(tf.concat([X_x_b1_train, X_y_b1_train], 1), weights, biases)

    # 边界条件预测 (y = ya, yb)
    u_y_b0_pred = pinn.fnn(tf.concat([Y_x_b0_train, Y_y_b0_train], 1), weights, biases)
    u_y_b1_pred = pinn.fnn(tf.concat([Y_x_b1_train, Y_y_b1_train], 1), weights, biases)

    # 定义损失函数
    loss = 0.1 * tf.reduce_mean(tf.square(f_pred)) + \
           tf.reduce_mean(tf.square(u_x_b0_pred)) + \
           tf.reduce_mean(tf.square(u_x_b1_pred)) + \
           tf.reduce_mean(tf.square(u_y_b0_pred)) + \
           tf.reduce_mean(tf.square(u_y_b1_pred))

    # 定义优化器
    train_adam = tf.train.AdamOptimizer(0.0008).minimize(loss)
    train_lbfgs = tf.contrib.opt.ScipyOptimizerInterface(loss,
                                                         method="L-BFGS-B",
                                                         options={'maxiter': 10000 ,
                                                                  'ftol': 1.0 * np.finfo(float).eps})

    # 创建 TensorFlow 会话
    session = tf.Session()
    session.run(tf.global_variables_initializer())
    
    # 创建训练字典
    train_feed_dict = {x_res_train: X_res[:, 0:1], y_res_train: X_res[:, 1:2], 
                       X_x_b0_train: X_b0[:, 0:1], X_y_b0_train: X_b0[:, 1:2],
                       X_x_b1_train: X_b1[:, 0:1], X_y_b1_train: X_b1[:, 1:2], 
                       Y_x_b0_train: Y_b0[:, 0:1], Y_y_b0_train: Y_b0[:, 1:2], 
                       Y_x_b1_train: Y_b1[:, 0:1], Y_y_b1_train: Y_b1[:, 1:2]}

    # 创建训练模型
    model = Train(train_feed_dict)
    
    # 记录训练时间
    start_time = time.perf_counter()
    model.nntrain(session, u_pred, loss, train_adam, train_lbfgs)
    stop_time = time.perf_counter()
    print('训练时间为 %.3f 秒' % (stop_time - start_time))

    # 保存预测数据和图像
    num_test_x_points = 101
    num_test_y_points = 101
    data_saver = SavePlot(session, x_range, y_range, num_test_x_points, num_test_y_points, xa, xb, ya, yb)
    data_saver.save_and_plot(u_pred, x_res_train, y_res_train)

    
        # 计算误差
    x_test = np.linspace(x_range[0], x_range[1], num_test_x_points).reshape((-1, 1))
    y_test = np.linspace(y_range[0], y_range[1], num_test_y_points).reshape((-1, 1))
    x_t, y_t = np.meshgrid(x_test, y_test)
    x_t = np.reshape(x_t, (-1, 1))
    y_t = np.reshape(y_t, (-1, 1))
    test_dict = {x_res_train: x_t, y_res_train: y_t}
    u_test_pred = session.run(u_pred, feed_dict=test_dict)  # 预测在均匀网格上的解
    Exact_sln = Exact(xa, xb, ya, yb)
    u_test_exact = Exact_sln.u_exact(np.hstack((x_t, y_t)))

    # 计算误差
    L2_error, max_error = compute_errors(u_test_pred, u_test_exact)
    print('L2 Error: %.6e' % L2_error)
    print('Max Error: %.7e' % max_error)
    
    
if __name__ == '__main__':
    main()
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129

程序中已经写好了详细的注释,关于优化器与 TF 会话(session) 的相关知识请各位移步 TensorFlow 优化器使用。另外建议读者对比阅读我之前总结的一维PINN 算法的实现 ,理解一维二维的本质不同,更高维的 PDE 求解也就不在话下了。

运行结果

在这里插入图片描述
在这里插入图片描述
效果不错!

-----------------------------------------------------------------------------------------------
本专栏目标从简单的一维 Poisson 方程,到对流扩散方程,Burges 方程,到二维,三维以及非线性方程,发展方程,积分方程等等,所有文章包含全部可运行代码。请持续关注!


作者 :计算小屋
个人主页计算小屋的主页

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/正经夜光杯/article/detail/889483
推荐阅读
相关标签
  

闽ICP备14008679号