python差分方程拟合 python求解差分方程

Python 一维波动方程数值解及可视化

一、效果展示

2. 两端固定，初值条件为 $\varphi(x) = \sin(3 \pi x)$



3. 右端为自由端，在前两秒施加外力，随后转为固定端

3. 两端施加不同频率外力

二、 求解原理

a. 微分方程

一维波动方程的一般形式如下
$\begin{cases} \dfrac{\partial^2u}{\partial t^2} = a^2\dfrac{\partial^2 u}{\partial x^2},0 <x <l, t > 0\\ \end{cases} \tag{1}$

b. 差分方程

我们先不考虑初值条件与边界条件，为了在不求该方程解析解的情况下描述方程图像，我们对原始方程进行差分处理。

设 $x, t$ 在数轴上被均匀的分割为 $n_x, n_t$等分段，每一段长度为$\Delta x, \Delta t$, 则第$i$段位移在第$j$ 段时间内，可以表示为$u_{i, j} = u(x_i, t_j)$

当 $\lim n_x \to \infty$ 且 $\lim n_t \to \infty$ 时，可以认为
$\frac{\partial u(x_{i}, t_j)}{\partial x} = \frac{u(x_{i+1}, t_j)-u(x_i, t_j)}{\Delta x} + o(\Delta x)\tag{2}$
进而有
$\begin{aligned} \frac{\partial^2{u}}{\partial x^2}(x_i, t_j) &= \frac{\dfrac{\partial u(x_{i+1}, t_j)}{\partial x} - \dfrac{\partial u(x_i, t_j)}{\partial x}}{\Delta x}\\ &=\frac{u_{i+1,j} -2u_{i, j} + u_{i-1, j}}{\Delta x^2} \end{aligned} \tag{3}$
同理可得
$\frac{\partial^2 u}{\partial t^2}(x_i, t_j) = \frac{u_{i, j + 1} - 2u_{i, j} + u_{i, j - 1}}{\Delta t^2}\tag{4}$
忽略高阶项，并将(3), (4)代入方程(1)中，得
$u_{i, j + 1} = (a\frac{\Delta t}{\Delta x})^2(u_{i + 1, j} -2u_{i, j} + u_{i - 1,j}+ 2u_{i, j} - u_{i, j - 1})$
上式中可以看到，如果需要求时间第$j + 1$时刻的波动方程，只需要知道$j, j - 1$ 时刻的$u$

我们将上式称为一维波动方程的差分形式。

c. 边界条件处理

1. 第一类边界条件

$u|_{x=0} = 0, u|_{l=0} = 0$

在差分方程中可以写成如下形式
$u_{1, j} = 0\\ u_{n_x,j} = 0$

2. 第二类边界条件
   $\left.\frac{\partial u}{\partial x}\right|_{x=0} = 0, \left.\frac{\partial u}{\partial x}\right|_{l=0} = 0,$
   差分形式:
   $\begin{aligned} u_{1, j}-u_{0, j} &= 0\\ u_{n_x, j} - u_{n_x-1, j} &= 0 \end{aligned}$
3. 第三类边界条件
   $\frac{\partial{u}}{\partial x} + \sigma u =f$
   差分形式:
   $\frac{u_{i, j} - u_{i - 1, j}}{\Delta x} + u_{i, j}= f$

三、Python 实现

1. 变量设计

# 波速
wave_velocity = 1
# x范围和节点数
x_min, x_max, number_dx = 0, 1, 100
# t范围和节点数
t_min, t_max, number_dt = 0, 2, 300
# x分片
x_ticks = np.linspace(x_min, x_max, number_dx)
t_ticks = np.linspace(t_min, t_max, number_dt)

# 每小段dx, dt的长度
dx = (x_max - x_min) / (number_dx - 1)
dt = (t_max - t_min) / (number_dt - 1)

# 时间步控制变量
time_step = 1

# 微分方程记录数组
u = np.zeros((number_dx, number_dt))

2. 初始化

phi = lambda x: np.sin(6 * np.pi * x)
for i, x in enumerate(self.x_ticks):
    # t = 0 时的初值
	u[i][0] = phi(x)

3. 迭代函数

def next_step():
    global time_step
    # 时间还没有结束时
    if t_ticks[self._time_step] < t_max:
        # 获取当前时间的波动方程各点值
        u_t = u[:, time_step]
        end = len(u_t) - 1
        # 计算下一个时间点的各点坐标值
        ddu = list(u_t[0: end - 1] - 2 * u_t[1: end] + u_t[2: end + 1])
        # 由于计算时，左右两个边界点无法通过差分计算(超过数值范围)
        # 这里将左右两个点先设置成0，之后再代入边界条件进行处理
        ddu = np.array([0] + ddu + [0])
        u[:, time_step + 1] = (wave_velocity * dt / dx) ** 2 * ddu + 2 * u[:, time_step] - u[:, time_step - 1]

        # 时间步自增，并且计算边界点的值
        # 这里左右端点都使用第一类边界条件
        time_step += 1
        u[0][time_step] = 0
		u[-1][time_step] = 0
        return time_step - 1
    else:
        # 此时迭代结束
        return -1

4. 执行

while True:
    result = next_step()
    if result == -1:
        break

5. 绘制

u_max = np.max(u)
for i in range(number_dt):
    plt.clf()
    plt.plot(x_ticks, u[:, i])
    plt.axis((x_min - x_max / 10, x_max + x_max / 10, -1.2 * u_max, 1.2 * u_max))
    plt.xlabel("Distance (x), t = {:.2f}(s)".format(t_ticks[i]))
    plt.ylabel("u")
    plt.pause(dt / 2)
plt.show()

6. 汇总

import numpy as np
import matplotlib.pyplot as plt


# 波速
wave_velocity = 1
# x范围和节点数
x_min, x_max, number_dx = 0, 1, 100
# t范围和节点数
t_min, t_max, number_dt = 0, 2, 300
# x分片
x_ticks = np.linspace(x_min, x_max, number_dx)
t_ticks = np.linspace(t_min, t_max, number_dt)

# 每小段dx, dt的长度
dx = (x_max - x_min) / (number_dx - 1)
dt = (t_max - t_min) / (number_dt - 1)

# 时间步控制变量
time_step = 1

# 微分方程记录数组
u = np.zeros((number_dx, number_dt))

phi = lambda x: np.sin(3 * np.pi * x)
for i, x in enumerate(x_ticks):
    # t = 0 时的初值
    u[i][0] = phi(x)


def next_step():
    global time_step 
    # 时间还没有结束时
    if t_ticks[time_step] < t_max:
        # 获取当前时间的波动方程各点值
        u_t = u[:, time_step]
        end = len(u_t) - 1
        # 计算下一个时间点的各点坐标值
        ddu = list(u_t[0: end - 1] - 2 * u_t[1: end] + u_t[2: end + 1])
        # 由于计算时，左右两个边界点无法通过差分计算(超过数值范围)
        # 这里将左右两个点先设置成0，之后再代入边界条件进行处理
        ddu = np.array([0] + ddu + [0])
        u[:, time_step + 1] = (wave_velocity * dt / dx) ** 2 * ddu + 2 * u[:, time_step] - u[:, time_step - 1]

        # 时间步自增，并且计算边界点的值
        # 这里左右端点都使用第一类边界条件
        time_step += 1
        u[0][time_step] = 0
        u[-1][time_step] = 0
        return time_step - 1
    else:
        # 此时迭代结束
        return -1


def main():
    while True:
        result = next_step()
        if result == -1:
            break

    u_max = np.max(u)
    for t in range(number_dt):
        plt.clf()
        plt.plot(x_ticks, u[:, t])
        plt.axis((x_min - x_max / 10, x_max + x_max / 10, -1.2 * u_max, 1.2 * u_max))
        plt.xlabel("Distance (x), t = {:.2f}(s)".format(t_ticks[t]))
        plt.ylabel("u")
        plt.pause(dt / 2)
    plt.show()


if __name__ == '__main__':
    main()

四、 存在的问题及改进

问题概述

1. 在上面的代码中，如果 number_dx 和 number_dt 设置不当，容易造成溢出。需要多次尝试设置这两个参数值。
2. 每次修改边界条件时较为麻烦。
3. 当number_dt很大时，生成的图像波动速读很慢，观感很差。
4. 如果需要在波动的过程中附加其他的外力等条件，需要另外修改代码。

改进方案

1. 自动的根据需求调整number_dx和number_dt，以防止溢出。
2. 将边界条件封装起来，传入所需要的第几类边界条件时可以自动更改。
3. 使用matplotlib 自带的动画生成函数，并封装帧数、延迟等的调整，使动画美观。
4. 提供callbacks 借口，在每次迭代后，把各点坐标的数组丢入函数中，用户可以设计callback函数实现丰富的功能。

五、改进后的实现

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as amt


class OneDimensionalFluctuation:
    def __init__(self,
                 *args,
                 length=1,
                 wave_velocity=1,
                 number_split=75,
                 end_time=2.495,
                 phi=lambda x: 0,
                 varphi=lambda x: 0,
                 callbacks=None,
                 left_boundary_situation=1,
                 right_boundary_situation=1,
                 ):
        self.wave_velocity = wave_velocity

        # init basic variables
        self.x_min, self.x_max, self.number_dx = 0, length, number_split  # Discretization of x
        self.t_min, self.t_max, self.number_dt = 0, end_time, int(number_split * end_time)  # Discretization of y

        self.x_ticks = np.linspace(self.x_min, self.x_max, self.number_dx)
        self.t_ticks = np.linspace(self.t_min, self.t_max, self.number_dt)

        self.dx = (self.x_max - self.x_min) / (self.number_dx - 1)
        self.dt = (self.t_max - self.t_min) / (self.number_dt - 1)

        self._time_step = 1

        # init functions
        self.phi = phi
        self.varphi = varphi
        self.callbacks = callbacks
        self.left_boundary_situation = left_boundary_situation
        self.right_boundary_situation = right_boundary_situation

        # the fluctuation
        self.u = np.zeros((self.number_dx, self.number_dt))

        # init u(x, 0)
        for i, x in enumerate(self.x_ticks):
            self.u[i][0] = phi(x)

    def _apply_boundary_situation(self):
        # left
        if self.left_boundary_situation == 1:
            self.u[0][self._time_step] = 0
        elif self.left_boundary_situation == 2:
            self.u[0][self._time_step] = self.u[1][self._time_step]
        else:
            raise ValueError("No such left boundary situation {}".format(self.left_boundary_situation))

        # right
        if self.right_boundary_situation == 1:
            self.u[-1][self._time_step] = 0
        elif self.right_boundary_situation == 2:
            self.u[-1][self._time_step] = self.u[-2][self._time_step]
        else:
            raise ValueError("No such right boundary situation {}".format(self.right_boundary_situation))

    def _next_step(self):
        if self.t_ticks[self._time_step] < self.t_max:
            # calculate each position's next value
            u_t = self.u[:, self._time_step]
            end = len(u_t) - 1
            ddu = list(u_t[0: end - 1] - 2 * u_t[1: end] + u_t[2: end + 1])
            ddu = np.array([0] + ddu + [0])
            self.u[:, self._time_step + 1] = (self.wave_velocity * self.dt / self.dx) ** 2 * ddu + 2 * self.u[:, self._time_step] - self.u[:, self._time_step - 1]

            # apply the boundary situation
            self._time_step += 1
            self._apply_boundary_situation()

            return self._time_step - 1
        else:
            # the iter is over
            return -1

    def run(self):
        while self._time_step < self.number_dt:
            if self.callbacks is not None:
                self.callbacks(self.t_ticks[self._time_step], self.dx, self.u[:, self._time_step])
            result = self._next_step()
            if result == -1:
                break

    def draw(self):
        # get the max u value
        u_max = np.max(self.u)
        for i in range(self.number_dt):
            plt.clf()
            plt.plot(self.x_ticks, self.u[:, i])
            plt.axis((self.x_min - self.x_max / 10, self.x_max + self.x_max / 10, -1.2 * u_max, 1.2 * u_max))
            plt.xlabel("Distance (x), t = {:.2f}(s)".format(self.t_ticks[i]))
            plt.ylabel("u")
            plt.pause(self.dt / 2)
        plt.show()

    def get_animation(self, name, fps=10, interval=None):
        if interval is None:
            interval = self.number_dt / (self.t_max + 1)

        fig, ax = plt.subplots()
        line, = ax.plot(self.x_ticks, self.u[:, 0])
        ax.set_xlim(self.x_min - self.x_max / 10, self.x_max + self.x_max / 10)
        ax.set_ylim(-1.1 * np.max(self.u), 1.1 * np.max(self.u))

        def update(num):
            line.set_ydata(self.u[:, num])
            ax.set_xlabel("t={:.2f}(s)".format(self.t_ticks[num]))

            return line,

        ani = amt.FuncAnimation(fig, update, frames=self.number_dt, interval=interval)
        ani.save('{}.gif'.format(name), fps=fps)
        plt.show()

        return ani

使用方法

1. 初始化参数列表

def __init__(self,
             *args,
             length=1,
             wave_velocity=1,
             number_split=75,
             end_time=2.495,
             phi=lambda x: 0,
             varphi=lambda x: 0,
             callbacks=None,
             left_boundary_situation=1,
             right_boundary_situation=1,
             )

参数名	含义	类型	默认值
length	弦的长度	float	1
wave_velocity	波速	float	1
number_split	节点的分割数(自动防溢出)	int	75
end_time	波动终止时刻	float	2.495
phi	位移初值条件	function	0
varphi	速度初值条件 (目前无效)	function	0
callbacks	回馈函数	function	None
left_boundary_situation	左边界条件	int	1
right_boundary_situation	右边界条件	int	1

1. 使用时，边界条件只能传入1 或 2，第三类边界条件还没有实现。
2. varphi 是速度初值条件，目前还没有实现。

2. 一些使用例子

1. 两端固定，初值条件为 $\varphi(x) = \sin(3 \pi x)$

if __name__ == '__main__':
    u = OneDimensionalFluctuation(
        length=1,
        wave_velocity=1,
        number_split=100,
        end_time=2,
        phi=lambda x: np.sin(3 * np.pi * x),
        left_boundary_situation=1,
        right_boundary_situation=1
        )
    u.run()
    u.draw()
    # u.get_animation("test")

1. 左端固定，右端前2s附加外力

def add_right_power(current_time, dx, ut):
    if current_time < 2:
        ut[-1] = ut[-2] + 10 * dx * np.sin(4 * np.pi * current_time)

        
if __name__ == '__main__':
   u = OneDimensionalFluctuation(
       callbacks=add_right_power, 
       end_time=np.pi
   )
    u.run()
    u.draw()
    # u.get_animation("right_power", fps=20, interval=70)

3. 右端为自由端，初值条件为 $\varphi(x) = \sin(6 \pi x)$

if __name__ == "__main__":
    phi = lambda x: np.sin(6 * np.pi * x)
 u = OneDimensionalFluctuation(phi=phi, end_time=3, length=1, right_boundary_situation=2, callbacks=add_right_power)
    u.run()
    u.draw()
    # u.get_animation("secondary boundary", fps=20, interval=100)

3. 左右两端为自由端，两端施加同频率等大外力

def add_two_power(current_time, dx, ut):
    if current_time < 2:
        ut[-1] = ut[-2] + 10 * dx * np.sin(4 * np.pi * current_time)
        ut[0] = ut[1] + 10 * dx * np.sin(4 * np.pi * current_time)


        
if __name__ == '__main__':
	u = OneDimensionalFluctuation(
        length=1,
        wave_velocity=1,
        number_split=100,
        end_time=2,
        callbacks=add_two_power,
        left_boundary_situation=1,
        right_boundary_situation=1,
        )
    u.run()
    u.draw()
    # u.get_animation("two_power")

5. 左右两端为自由端，两端施加不同频率外力

def add_two_power(current_time, dx, ut):
    if current_time < 2:
        ut[-1] = ut[-2] + 10 * dx * np.sin(3 * np.pi * current_time)
        ut[0] = ut[1] + 10 * dx * np.sin(4 * np.pi * current_time)


if __name__ == '__main__':
   u = OneDimensionalFluctuation(
        length=1,
        wave_velocity=1,
        number_split=100,
        end_time=1.9,
        callbacks=add_two_power,
    )
    u.run()
    u.draw()

附录

OneDimensionalFluctuation API

Attributes

1. wave_velocityint the speed of the wave.
2. x_min, x_max, number_dxfloat, float, int the minimum, maximum of x and the segment number of x.
3. t_min, t_max, number_dt
   float, float, int the minimum, maximum of t and the segment number of t.
4. x_ticks
   List[float] the value of x in each segment.
5. t_ticks
   List[float] the value of t in each segment.
6. dx
   float the distance between each x segment
7. dt
   float the distance between each t segment.
8. _time_step
   int[protected] the recorder of time step.
9. uList[List[float]] the array that save the value of the wave.
   u[i][j] represent the value of the point at the i-th segment in x_ticks when the time is j-th of t_ticks.
10. phi
    function(float) -> float this function is the initial function of u(x, 0), will give the value of x position, and need to return the value of u(0, x).
11. varphi
    function(float) -> float this function is the initial function of x, will give the value of x position, and need to return the value of u(0, x).
    this attribute is no use now
12. callbacks
    function(float, float, List[float]) After each iteration, the callbacks function will be called, and each of the argument are: current_iteration_time, dx, ut.
    ut is a array of the wave at current iteration, and you can change this array to change the wave.

Example:

def add_right_power(current_time, dx, ut):
    if current_time < 2:
        ut[-1] = ut[-2] + 10 * dx * np.sin(4 * np.pi * current_time)

This function will add a power to the right point of the wave if the current time is less than 2(s).

1. left_boundary_situation
   int only two option: 1 and 2.
   1: first boundary situation.
   2: second boundary situation.
2. right_boundary_situation
   int only two option: 1 and 2.
   1: first boundary situation.
   2: second boundary situation.

Method

1. run
   solve the problem, you need to run this function before draw and get animation.
2. draw
   draw a gif by matplotlib.pyplot.plc().
3. get_animation
   show an animation by matplotlib.animation. And will return the Animation Object.

参考

1. 一维波动方程的数值解