pinn求解常微分方程

### 使用PINN求解常微分方程

对于常微分方程（ODE），物理信息神经网络（PINNs）提供了一种新颖而强大的数值求解方法。这种方法的核心在于构建一个能够满足给定ODE及其边界条件的神经网络模型。

为了使神经网络逼近ODE的解，必须设计合适的损失函数来反映这些约束。具体来说，损失函数由两部分组成：一部分基于ODE本身，另一部分则考虑初始或边界条件[^1]。

#### 构建神经网络架构

通常情况下，会选择一个多层感知机作为基础结构。输入节点接收自变量（例如时间t），输出则是待求解函数u(t)。隐藏层数量及每层中的神经元数目依据问题复杂度调整。

import tensorflow as tf
from tensorflow.keras import layers, models

def create_model(input_dim=1, output_dim=1):
    model = models.Sequential()
    model.add(layers.InputLayer(input_shape=(input_dim,)))
    
    # 添加多个全连接层
    for units in [20]*8:  # 这里设置为8层，每层有20个单元
        model.add(layers.Dense(units, activation='tanh'))
        
    model.add(layers.Dense(output_dim))
    return model

#### 定义残差项与损失函数

定义好网络之后，下一步就是计算ODE对应的残差项，并将其加入到总损失之中。这里假设有一个简单的二阶线性齐次ODE：

\[ u''(t)+\omega^{2}u(t)=f(t), \quad t\in[a,b], \]

其中\( f(t)\equiv 0\)表示无源项情况下的自由振动；另外还需要指定两个端点处的具体取值或者导数关系作为辅助条件。

通过自动微分技术可以方便地获取任意阶导数表达式，在TensorFlow框架下这一步骤变得尤为简单。

@tf.function
def compute_residuals(model, t, omega_squared):
    with tf.GradientTape(persistent=True) as tape:
        tape.watch(t)
        u_pred = model(t)

        du_dt = tape.gradient(u_pred, t)
    ddu_ddt = tape.gradient(du_dt, t)
    
    residuals = ddu_ddt + omega_squared * u_pred
    
    del tape
    return residuals

最后组合上述各项形成完整的损失公式并利用优化算法最小化之以获得最优参数配置。

optimizer = tf.optimizers.Adam()

for epoch in range(num_epochs):
    with tf.GradientTape() as tape:
        pred_u = model(train_t)
        resids = compute_residuals(model, train_t, omega_sqrd)
        
        bc_loss = ...  # 边界/初值条件造成的误差平方和
        
        pde_loss = tf.reduce_mean(tf.square(resids))  
        total_loss = pde_loss + bc_loss
        
    grads = tape.gradient(total_loss, model.trainable_variables)
    optimizer.apply_gradients(zip(grads, model.trainable_variables))

    if (epoch % display_step == 0): print(f'Epoch {epoch}, Loss={total_loss.numpy()}')