PDE保能量数值格式有什么意义

最新推荐文章于 2024-07-24 22:13:06 发布

懒鬼|

最新推荐文章于 2024-07-24 22:13:06 发布

阅读量984

点赞数 26

文章标签：算法数学建模

本文链接：https://blog.csdn.net/m0_60307945/article/details/139720688

版权

在物理学中，"保能量"通常指的是一个系统在没有外部能量输入或损失的情况下，其总能量保持不变。例如，在经典力学中，如果一个系统是封闭的（没有能量的流入或流出），那么这个系统的总能量（动能加势能）将保持常数。对于方程的解，如果一个方程或系统保持某种物理量（如能量、动量等）不变，这通常会对方程的解的性质有重要影响。比如在物理中，保能量的解法可能会确保数值解的长期行为（如稳定性和误差控制）更加符合实际物理情况。

在数学中，对于同一个方程组系统，使用保能量的数值解法和不保能量的数值解法得到的结果确实可能会有所不同，特别是在长时间模拟中更为明显。这些差异主要体现在以下几个方面：

长期行为： 保能量解法在长时间模拟中能够更好地维持系统的物理或数学性质，如能量守恒。这对于需要长时间观察系统行为的研究特别重要，比如在天体动力学和量子力学模拟中。不保能量的方法可能导致能量逐渐累积或耗散，从而引入长期的计算误差。
数值稳定性： 保能量方法往往能提供更好的数值稳定性，尤其是在处理保守系统（系统的总能量应保持不变的系统）时。不保能量的数值方法可能在一定时间后由于积累的数值误差导致解发散或偏离真实解。
物理合理性： 在某些物理问题中，确保数值解法能保持能量守恒是非常重要的，以确保数值解的物理合理性。而不保能量的解法可能在模拟结果中引入非物理的行为，这在解释或利用这些结果时可能带来问题。
解的精度和质量： 保能量的解法通常在维持解的精度和质量方面表现更好，尤其是在需要精确遵循物理定律的模拟中。不保能量的方法可能在较短的时间内给出接近的结果，但随着时间的推移，误差可能逐渐显现。
特定应用的要求： 根据具体应用的需要，选择保能量还是不保能量的数值方法也有所不同。例如，在需要极高精度的科学研究或工程应用中，可能优先选择保能量方法。

接下来考虑一个简单的物理问题：一个理想的单摆系统。单摆系统是一个非常经典的例子，可以很好地展示保能量和不保能量数值解法的效果差异。

问题描述

理想单摆由一个无质量的绳子和一个质点组成，绳子的一端固定，质点在重力作用下摆动。其运动方程可以通过牛顿第二定律导出，也可以通过哈密顿动力学得到。在理想单摆的情况下，系统的总能量（动能加势能）应当是守恒的，即：

$E= \frac{1}{2}mL^2\omega^2+mgL(1-\cos(\theta))$ .

数值方法如果能保持这个能量不变，那么模拟结果通常更为可靠，特别是在需要长时间观察系统行为时。辛积分方法通过保持系统哈密顿结构的方式来实现能量守恒，而常规的欧拉方法则可能因为数值误差的累积而导致能量逐渐偏离其初始值。这里我们使用角度 $\theta$ 作为描述状态的变量，单摆的运动方程可以表示为：

$\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2}+\frac{g}{L}\sin(\theta)=0$ ,

其中，g 是重力加速度，L 是摆绳的长度。定义角速度 $\omega=\frac{\mathrm{d}\theta}{\mathrm{d}t}$ , 则可以将上述二阶常微分方程转化为两个一阶方程，即：

$\begin{cases} \frac{\mathrm{d}\theta}{\mathrm{d}t}=\omega, \\ \frac{\mathrm{d}\omega}{\mathrm{d}t}=-\frac{g}{L}\sin(\theta). \end{cases}$

我们将使用两种方法来解这个方程：

欧拉方法（不保能量）： 这是最简单的数值积分方法之一，但它不保证能量守恒。对于上述的两个方程，更新规则为：

$\theta_{n+1}=\theta_n+\omega_n \cdot \mathrm{d}t, \\ \omega_{n+1}=\omega_n-\frac{g}{L}\sin(\theta_n) \cdot \mathrm{d}t.$
这种方法是显式的，容易实现，但它不保守能量，尤其是在长时间积分和较大的时间步长下，能量的计算误差会累积。
辛欧拉方法（保能量）： 辛欧拉方法是辛积分方法的一种，它特别适用于解哈密顿系统。辛方法的特点是在数值积分过程中保持系统的哈密顿结构，从而保守能量。对于单摆问题，辛欧拉方法的更新规则可以表述为：

$\omega_{n+1}=\omega_n-\frac{g}{L}\sin(\theta_n) \cdot \mathrm{d}t, \\ \theta_{n+1}=\theta_n+\omega_{n+1} \cdot \mathrm{d}t.$
这种方法先更新角速度，然后用新的角速度来更新角度，这有助于更好地保守能量。

实现

我们将模拟单摆在没有阻尼和外部力的条件下的运动，初始时刻摆偏离垂直位置一个小角度，比如 5°，没有初始角速度。

以下是Python代码实现的简要说明，包括两种方法的比较：

import numpy as np
import matplotlib.pyplot as plt

# 定义常数
g = 9.81  # 重力加速度, m/s^2
L = 1.0   # 摆长, m
m = 1.0   # 质量, kg (假设为1kg以简化计算)
dt = 0.01 # 时间步长, s
T = 10    # 总模拟时间, s

# 初始条件
theta_0 = np.deg2rad(5)  # 初始角度转换为弧度
omega_0 = 0  # 初始角速度

# 时间数组
t = np.arange(0, T+dt, dt)

# 初始化数组
theta_euler = np.zeros_like(t)
omega_euler = np.zeros_like(t)
energy_euler = np.zeros_like(t)
theta_euler[0] = theta_0
omega_euler[0] = omega_0

theta_symplectic = np.zeros_like(t)
omega_symplectic = np.zeros_like(t)
energy_symplectic = np.zeros_like(t)
theta_symplectic[0] = theta_0
omega_symplectic[0] = omega_0

# 欧拉方法积分
for i in range(1, len(t)):
    omega_euler[i] = omega_euler[i-1] - (g/L) * np.sin(theta_euler[i-1]) * dt
    theta_euler[i] = theta_euler[i-1] + omega_euler[i-1] * dt
    energy_euler[i] = 0.5 * m * L**2 * omega_euler[i]**2 + m * g * L * (1 - np.cos(theta_euler[i]))

# 辛欧拉方法积分
for i in range(1, len(t)):
    omega_symplectic[i] = omega_symplectic[i-1] - (g/L) * np.sin(theta_symplectic[i-1]) * dt
    theta_symplectic[i] = theta_symplectic[i-1] + omega_symplectic[i] * dt
    energy_symplectic[i] = 0.5 * m * L**2 * omega_symplectic[i]**2 + m * g * L * (1 - np.cos(theta_symplectic[i]))

# 绘制角度变化
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t, theta_euler, label='Euler Method')
plt.plot(t, theta_symplectic, label='Symplectic Euler Method')
plt.xlabel('Time (s)')
plt.ylabel('Theta (rad)')
plt.title('Angle Over Time')
plt.legend()

# 绘制能量变化
plt.subplot(1, 2, 2)
plt.plot(t, energy_euler, label='Euler Method')
plt.plot(t, energy_symplectic, label='Symplectic Euler Method')
plt.xlabel('Time (s)')
plt.ylabel('Energy (Joules)')
plt.title('Energy Over Time')
plt.legend()

plt.tight_layout()
plt.show()

结果如下图所示，欧拉方法的能量逐渐变大，辛方法的能量虽然有轻微波动，但总体保持不变。

以下，我还尝试了使用向后欧拉方法来求解，附代码及结果，此时向后欧拉方法的能量递减。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# 定义常数
g = 9.81  # 重力加速度, m/s^2
L = 1.0   # 摆长, m
m = 1.0   # 质量, kg
dt = 0.01 # 时间步长, s
T = 10    # 总模拟时间, s

# 初始条件
theta_0 = np.deg2rad(5)  # 初始角度转换为弧度
omega_0 = 0  # 初始角速度

# 时间数组
t = np.arange(0, T+dt, dt)

# 初始化数组
theta_back = np.zeros_like(t)
omega_back = np.zeros_like(t)
energy_back = np.zeros_like(t)
theta_back[0] = theta_0
omega_back[0] = omega_0

theta_sym = np.zeros_like(t)
omega_sym = np.zeros_like(t)
energy_sym = np.zeros_like(t)
theta_sym[0] = theta_0
omega_sym[0] = omega_0

# 向后欧拉方法积分
def backward_euler(y, y_prev):
    theta_prev, omega_prev = y_prev
    theta, omega = y
    return [theta - theta_prev - omega * dt, omega - omega_prev + (g/L) * np.sin(theta) * dt]

for i in range(1, len(t)):
    y_prev = [theta_back[i-1], omega_back[i-1]]
    y = fsolve(backward_euler, y_prev, args=(y_prev,))
    theta_back[i], omega_back[i] = y
    energy_back[i] = 0.5 * m * L**2 * omega_back[i]**2 + m * g * L * (1 - np.cos(theta_back[i]))

# 辛欧拉方法积分
for i in range(1, len(t)):
    omega_sym[i] = omega_sym[i-1] - (g/L) * np.sin(theta_sym[i-1]) * dt
    theta_sym[i] = theta_sym[i-1] + omega_sym[i] * dt
    energy_sym[i] = 0.5 * m * L**2 * omega_sym[i]**2 + m * g * L * (1 - np.cos(theta_sym[i]))

# 绘制结果
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(t, theta_back, label='Backward Euler Method')
plt.plot(t, theta_sym, label='Symplectic Euler Method', linestyle='--')
plt.xlabel('Time (s)')
plt.ylabel('Theta (rad)')
plt.title('Angle Over Time')
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(t, energy_back, label='Backward Euler Energy')
plt.plot(t, energy_sym, label='Symplectic Euler Energy', linestyle='--')
plt.xlabel('Time (s)')
plt.ylabel('Energy (Joules)')
plt.title('Energy Over Time')
plt.legend()

plt.tight_layout()
plt.show()

结论

系统的能量守恒或耗散

系统是否具有能量守恒或耗散的特性，是由系统的物理或数学本质决定的。例如：

在保守系统中，如完美的单摆系统，没有外部力和摩擦的情况下，系统的总能量（动能和势能之和）应该是守恒的。
在非保守系统中，如有阻尼的振子，系统能量随时间逐渐减少，表现为能量的耗散。

这些特性是由系统的基本方程（如牛顿第二定律、哈密顿系统方程等）决定的。

数值离散格式的性态保持

当我们使用数值方法来求解这些系统的方程时，理论上希望数值解能够尽可能地保留这些物理或数学特性。然而，由于数值方法的局限性（如截断误差、舍入误差等），往往数值解并不能完全保持系统的原始性态，如能量守恒。因此：

有些数值方法，如辛方法在求解哈密顿系统时，被特别设计以保持系统的哈密顿结构，从而尽可能保持能量守恒。
而像普通的欧拉方法、龙格-库塔方法等，尽管它们广泛应用于各种动力学系统的数值求解，但并不保证能量守恒，尤其是在长时间积分中。

数值方法的正确性证明

在开发或选择特定的数值离散格式时，通常需要进行理论分析或实验验证来确认这些方法是否能保持系统的原有性态：

理论分析：通过数学证明来分析数值方法是否保持如能量守恒、动量守恒等物理定律。例如，通过证明数值方法保持辛结构，可以推断其在哈密顿系统中的表现。
实验验证：通过数值实验，如模拟已知解的问题并与理论解进行比较，来观察数值解的行为是否符合预期，如能量是否守恒。

懒鬼|

关注

26
点赞
踩
29

收藏

觉得还不错? 一键收藏
0
评论
PDE保能量数值格式有什么意义

在物理学中，"保能量"通常指的是一个系统在没有外部能量输入或损失的情况下，其总能量保持不变。例如，在经典力学中，如果一个系统是封闭的（没有能量的流入或流出），那么这个系统的总能量（动能加势能）将保持常数。对于方程的解，如果一个方程或系统保持某种物理量（如能量、动量等）不变，这通常会对方程的解的性质有重要影响。比如在物理中，保能量的解法可能会确保数值解的长期行为（如稳定性和误差控制）更加符合实际物理情况。
复制链接

扫一扫