精华内容
下载资源
问答
  • Python 一维波动方程数值解及可视化 一、效果展示 两端固定,初值条件为 φ(x)=sin⁡(3πx)\varphi(x) = \sin(3 \pi x)φ(x)=sin(3πx) 右端为自由端,初值条件为 φ(x)=sin⁡(6πx)\varphi(x) = \sin(6 \pi x)...

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

    一、效果展示

    1. 两端固定,初值条件为 φ(x)=sin(3πx)\varphi(x) = \sin(3 \pi x)
      在这里插入图片描述

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

    在这里插入图片描述

    1. 两端施加不同频率外力
      在这里插入图片描述

    二、 求解原理

    a. 微分方程

    一维波动方程的一般形式如下
    {2ut2=a22ux2,0<x<l,t>0(1) \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,tx, t 在数轴上被均匀的分割为 nx,ntn_x, n_t等分段,每一段长度为Δx,Δt\Delta x, \Delta t, 则第ii段位移在第jj 段时间内,可以表示为ui,j=u(xi,tj)u_{i, j} = u(x_i, t_j)

    limnx\lim n_x \to \inftylimnt\lim n_t \to \infty 时,可以认为
    u(xi,tj)x=u(xi+1,tj)u(xi,tj)Δx+o(Δx)(2) \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}
    进而有
    2ux2(xi,tj)=u(xi+1,tj)xu(xi,tj)xΔx=ui+1,j2ui,j+ui1,jΔx2(3) \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}
    同理可得
    2ut2(xi,tj)=ui,j+12ui,j+ui1,jΔt2(4) \frac{\partial^2 u}{\partial t^2}(x_i, t_j) = \frac{u_{i, j + 1} - 2u_{i, j} + u_{i - 1, j}}{\Delta t^2}\tag{4}
    忽略高阶项,并将(3), (4)代入方程(1)中,得
    ui,j+1=(aΔtΔx)2(ui+1,j2ui,j+ui1,j+2ui,jui,j1) 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+1j + 1时刻的波动方程,只需要知道j,j1j, j - 1 时刻的uu 的函数值,因此,只要给出初值条件,该方程就可以通过差分递归求解。

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

    c. 边界条件处理

    1. 第一类边界条件

    ux=0=0,ul=0=0 u|_{x=0} = 0, u|_{l=0} = 0

    ​ 在差分方程中可以写成如下形式
    u1,j=0ul,j=0 u_{1, j} = 0\\ u_{l,j} = 0

    1. 第二类边界条件
      uxx=0=0,uxl=0=0, \left.\frac{\partial u}{\partial x}\right|_{x=0} = 0, \left.\frac{\partial u}{\partial x}\right|_{l=0} = 0,

      差分形式:
      u1,ju0,j=0ul,jul1,j=0 \begin{aligned} u_{1, j}-u_{0, j} &= 0\\ u_{l, j} - u_{l-1, j} &= 0 \end{aligned}

    2. 第三类边界条件
      ux+σu=f \frac{\partial{u}}{\partial x} + \sigma u =f
      差分形式:
      ui,jui1,jΔx+ui,j=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 时的初值
    	self.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_dxnumber_dt 设置不当,容易造成溢出。需要多次尝试设置这两个参数值。
    2. 每次修改边界条件时较为麻烦。
    3. number_dt很大时,生成的图像波动速读很慢,观感很差。
    4. 如果需要在波动的过程中附加其他的外力等条件,需要另外修改代码。

    改进方案

    1. 自动的根据需求调整number_dxnumber_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. 使用时,边界条件只能传入12,第三类边界条件还没有实现。

    2. varphi 是速度初值条件,目前还没有实现。

    2. 一些使用例子

    1. 两端固定,初值条件为 φ(x)=sin(3πx)\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")
      

      在这里插入图片描述

    2. 左端固定,右端前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)
      

    在这里插入图片描述

    1. 右端为自由端,初值条件为 φ(x)=sin(6πx)\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)
      

      在这里插入图片描述

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

      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")
      

    在这里插入图片描述

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

      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_velocity
      int
      the speed of the wave.

    2. x_min, x_max, number_dx
      float, 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. u
      List[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).

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

    14. 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. 一维波动方程的数值解
    展开全文
  • 色谱方程的有限差分求解器 该项目提供了功能和使用输运弥散模型模拟多个组件色谱的示例,并支持Hery,Langmuir和SMA作为组件的等温线模型。 使用一维有限差分Euler foward方法,使用numpy进行更快的计算,使用...
  • 关于斐波那契数列(Fibonacci sequence)之差分方程求解算法一点点记录! 斐波那契数列,又称之为黄金分割数列。具体指的是这样的一个数列: ​ 利用小学数学就能演算出后面的项数,可以很清楚的得到数列的后一项的...

    关于斐波那契数列(Fibonacci sequence)之差分方程求解算法一点点记录!

    斐波那契数列,又称之为黄金分割数列。具体指的是这样的一个数列:
    在这里插入图片描述
    利用小学数学就能演算出后面的项数,可以很清楚的得到数列的后一项的数值可以由其前面两项的到的。不过得注意的是最初的两项是第0个和第1个数值是分别是0和1,这样不难想到,肯定必须确定前面2项才能确定第三项的值,才能依次确定后面几项,具体规律可以用下面的这个通项式子表示

    在这里插入图片描述

    如果直接用递归的思想,估计很快就能写出来这个计数方法。这里我选用python做,工具嘛,哪个方便用哪个。

    def Fibonacci(N):
        if(N<2):
            return N
        else:
            return Fibonacci(N-1)+ Fibonacci(N-2)
    
    if __name__ == '__main__':
        with open("Fibonacci.txt","w") as f:
            for i in range(20):
                print(Fibonacci(i),end=" ")
                f.write(str(Fibonacci(i))+"  ")
            f.close()
    

    python Fibonacci.py 之后就会产生一个文件Fibonacci.txt里面内容如下:

    在这里插入图片描述

    这样迭代对于这样20数据来讲本来就没有压力。我简单测试一哈关于20,25,30,35,40次迭代完成,我电脑需要多少时间。

    在这里插入图片描述

    要是迭代100次的话,电脑不疯,我都要疯了,经过测试迭代40次左右的时间,我的台式电脑大致需要1分钟15秒左右。实在不敢测试100次的结果。
    在这里插入图片描述

    def Fibonacci(N):
        if(N<2):
            return N
        else:
            return Fibonacci(N-1)+ Fibonacci(N-2)
    value = [20,25,30,35,40]
    import time
    if __name__ == '__main__':
        for count in value:
            s = time.time()
            for i in range(count):
                print(Fibonacci(i),end=" ")
            print("\r\n花费时间为: ",(time.time() - s), "s")
    

    如何提高这个运算效率呢???

    解决的办法还是得落在这个表达式里面

    在这里插入图片描述

    这个形式就是差分方程的形式,我记得这是考研的数三的考点,奈何我考的是数一,不过学过微分方程,大致看了一哈,带公式计算即可,没有太复杂的东西。首先这个表达式是二阶常系数齐次线性差分方程。所谓二阶指的是该差分方程的最高阶数有表达式可以得出,所谓常系数指的是F(N+t) (t为任意自然数)前面的系数为常数,齐次指的是该方程右边为0。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    在这里插入图片描述

    斐波那契数列之差分方程求解过程

    在这里插入图片描述

    利用这个最后的表达式很快敲出代码:

    def Fibonacci(N):
        if(N<2):
            return N
        else:
            return Fibonacci(N-1)+ Fibonacci(N-2)
    
    def Fibonacci2(N):
        sqrt5 = 5**0.5
        fibnN = ((1+sqrt5)/2)**N-((1-sqrt5)/2)**N
        return round(fibnN/sqrt5)
    
    value = [20,25,30,35,40]
    import time
    
    if __name__ == '__main__':
        for count in value:
            s = time.time()
            for i in range(count):
                print(Fibonacci2(i),end=" ")
            print("\r\n花费时间为: ",(time.time() - s), "s")
      
    
    """
    使用Fibonacci函数
    输出:
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 
    花费时间为:  0.004997968673706055 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 
    花费时间为:  0.056983232498168945 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 
    花费时间为:  0.6223404407501221 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 1346269 2178309 3524578 5702887 
    花费时间为:  6.873619079589844 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 1346269 2178309 3524578 5702887 9227465 14930352 24157817 39088169 63245986 
    花费时间为:  75.94436860084534 s
    """
    
    
    """
    使用Fibonacci2函数
    输出:
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181
    花费时间为:  0.0 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 
    花费时间为:  0.0 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 
    花费时间为:  0.0 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 1346269 2178309 3524578 5702887 
    花费时间为:  0.0 s
    0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 1346269 2178309 3524578 5702887 9227465 14930352 24157817 39088169 63245986 
    花费时间为:  0.00099945068359375 s
    
    """
    

    尝试计算迭代100次,发现Fibonacci2算法仍旧快速输出结果。学好高数,真是能快速解决问题。不知道迭代次数足够多,是否有错误,没有验证过,但是计算是真的快啊

    每天学习亿点点,每天进步亿点点!(仅仅是记录自己的学习)要是有错误、或者其他简洁记得请指出来哈!

    展开全文
  • 一、题目 二、代码 from scipy.integrate import quad import numpy as np # 函数f(u)=1/2*u**2,故f`(u)=u def f_positive(upp_value): # 积分f+(u)中需要使用的函数 is_over_zero = int(upp_value >...

    一、题目
    1-1
    二、代码

    from scipy.integrate import quad
    import numpy as np
    
    
    # 函数f(u)=1/2*u**2,故f`(u)=u
    def f_positive(upp_value):
        # 积分f+(u)中需要使用的函数
        is_over_zero = int(upp_value > 0)
        return is_over_zero * upp_value
    
    
    def f_negative(upp_value):
        # 积分f-(u)中需要使用的函数
        is_over_zero = int(upp_value > 0)
        return (1 - is_over_zero) * upp_value
    
    
    def cal_fu(upp_value, form='+'):
        """
        :param upp_value:积分上限
        :param form: 可选值列表[‘+’,‘-’],决定返回f+(u)还是f-(u)
        :return: quad_r: f+(u)、f-(u)的值
        """
        quad_r = 0
        if form == "+":
            quad_r = quad(f_positive, 0, upp_value)[0]
        elif form == '-':
            quad_r = quad(f_negative, 0, upp_value)[0]
        return quad_r
    
    
    def cal_next_step(pt, px):
        """
        该空间点在下一时间层的值
        :param pt: 当前点在网格的空间位置
        :param px: 当前点在网格的时间位置
        :param grid: 网格点上的值
        :return:
        """
        u_j_n = grid_value[pt][px]  # 计算U(j,n)
        u_jplus1_n = grid_value[pt][px+1]  # 计算U(j+1,n)
        u_jminus1_n = grid_value[pt][px-1]  # 计算U(j-1,n)
        temp_minus = cal_fu(u_jplus1_n, form='-') - cal_fu(u_j_n, form='-')
        temp_plus = cal_fu(u_j_n) + cal_fu(u_jminus1_n)
        u_j_nplus1 = u_j_n - grid_ratio * (temp_minus + temp_plus)
        grid_value[pt+1][px] = u_j_nplus1
    
    
    if __name__ == "__main__":
        x_range = [-2, 2]  # 空间范围
        t_range = [0, 0.9]  # 时间范围
        delta_x = 0.1  # 空间步长
        delta_t = 0.01  # 时间步长
        grid_ratio = delta_t / delta_x  # 网格比
        grid_x = int((x_range[1] - x_range[0]) / delta_x) + 1  # 空间网格点数,此例中为41
        grid_t = int((t_range[1] - t_range[0]) / delta_t) + 1  # 时间网格点数,此例中为91
    
        # 考虑用列表grid_value来存储Ujn[[t=0.01],...,[t=0.9]]
        grid_value = np.zeros((grid_t, grid_x))  # 行代表某个时间、列代表某个空间
    
        # 将初始值t=0添加到grid_value中,即初始条件
        for i in range(grid_x):
            x_current = x_range[0] + delta_x * i
            if x_current > 0:
                grid_value[0][i] = 1
            else:
                grid_value[0][i] = -1
    
        # 将每一个时间层上的左右边界赋固定值
        grid_value[:, grid_x-1] = 1  # 右边界为1
        grid_value[:, 0] = -1  # 左边界为-1
    
        # 开始计算,时间上索引从0算到89,空间上索引从1算到39
        # 假设右边界必定收敛
        for i in range(0, grid_t-1):
            for j in range(1, grid_x-1):
                cal_next_step(i, j)
            grid_value[i, -1] = grid_value[i, -2]
        grid_value[-1, -1] = grid_value[-1, -2]
    
        # 仅将最后一个时间层的网格点数据保存到"2.txt"中
        np.savetxt('2.txt', grid_value[-1, :], fmt='%0.8f')
    

    三、运行结果
    3-1
    绘图如下:
    3-2

    展开全文
  • 这是用Python编写的一维模型,该模型使用有限差分求解半导体泊松漂移扩散方程。 该模型模拟了光照下的太阳能电池,但也可以适用于其他半导体器件。 可以对其进行修改以解决其他系统(即,通过更改边界条件,添加重组...
  • 该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组...

    有限差分法
    有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。
    该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。

    分类
    对于有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑,可分为中心格式和逆风格式。考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。差分方法主要适用于有结构网格,网格的步长一般根据实际情况和条件来决定。

    构造差分的方法
    构造差分的方法有多种形式,目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等,其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合,可以组合成不同的差分计算格式。

    泰勒级数
    在这里插入图片描述
    在这里插入图片描述

    Python中调用sympy模块求解差分

    from sympy import diff
    from sympy import symbols
    import sympy
    def func(t):
        return 2000 * sympy.log(14*10000/(14*10000-2100*t))-9.8*t  #函数
    t = symbols("t")
    print(func(16))
    print(diff(func(t),t))  #对函数进行求导
    print(diff(func(t),t).subs(t,16))
    

    运行上述程序结果:

    392.073691403521
    588000000000*(1 - 3*t/200)/(140000 - 2100*t)**2 - 9.8
    29.6736842105263
    

    【注意】:
    函数的写法不能调用sympy以外的模块的函数
    比如这样写就不对:

    from sympy import diff
    from sympy import symbols
    import math
    def func(t):
        return 2000 * math.log(14*10000/(14*10000-2100*t))-9.8*t
    t = symbols("t")
    print(func(16))
    print(diff(func(t),t))
    print(diff(func(t),t).subs(t,16))
    

    结果:

    392.07369140352074
    print(diff(func(t),t))
    return 2000 * math.log(14*10000/(14*10000-2100*t))-9.8*t
    raise TypeError("can't convert expression to float")
    TypeError: can't convert expression to float
    

    错误点在于:sympy的模块的函数不能调用math模块的函数

    一阶向前差分

    在这里插入图片描述
    Python代码:

    import sympy
    from sympy import diff
    from sympy import symbols
    
    #差分的对象
    x = 16
    k = 2  #步长
    x1 = x + k  #向前
    x2 = x - k  #向后
    
    #方程式
    def func(t):
        return 2000 * sympy.log(14*10000/(14*10000-2100*t))-9.8*t
    
    #一阶向前差分
    def for_difference():
        a_for_diff = (func(x1) - func(x))/k
        for_error = abs(a_for_diff - a_true)/a_true
        print(f'{x}的一阶向前差分值:{a_for_diff}')
        print(f'{x}的一阶向前差分的误差:{for_error*100}%')
    
    if __name__ == '__main__':
        t = symbols("t")
        a_true = diff(func(t), t).subs(t, x)  # 真值
        for_difference()
    

    结果:

    16的一阶向前差分值:30.4738991379398
    16的一阶向前差分的误差:2.69671578943888%
    

    一阶向后差分

    在这里插入图片描述
    Python代码:

    import sympy
    from sympy import diff
    from sympy import symbols
    
    #差分的对象
    x = 16
    k = 2  #步长
    x1 = x + k  #向前
    x2 = x - k  #向后
    
    #方程式
    def func(t):
        return 2000 * sympy.log(14*10000/(14*10000-2100*t))-9.8*t
    
    #一阶向后差分
    def beh_difference():
        a_beh_diff = (func(x) - func(x2))/k
        beh_error = abs(a_beh_diff - a_true)/a_true
        print(f'{x}的一阶向后差分值:{a_beh_diff}')
        print(f'{x}的一阶向后差分的误差:{beh_error*100}%')
    
    if __name__ == '__main__':
        t = symbols("t")
        a_true = diff(func(t), t).subs(t, x)  # 真值
        beh_difference()
    

    运行结果:

    16的一阶向后差分值:28.9145121806904
    16的一阶向后差分的误差:2.55840166138381%
    

    一阶中心差分

    在这里插入图片描述
    Python代码:

    import sympy
    from sympy import diff
    from sympy import symbols
    
    #差分的对象
    x = 16
    k = 2  #步长
    x1 = x + k  #向前
    x2 = x - k  #向后
    
    #方程式
    def func(t):
        return 2000 * sympy.log(14*10000/(14*10000-2100*t))-9.8*t
    
    #一阶中心差分
    def cen_difference():
        a_cen_diff = (func(x1)-func(x2))/(k*2)
        cen_error = abs(a_cen_diff - a_true)/a_true
        print(f'{x}的二阶中心差分值:{a_cen_diff}')
        print(f'{x}的二阶中心差分的误差:{cen_error*100}%')
    
    if __name__ == '__main__':
        t = symbols("t")
        a_true = diff(func(t), t).subs(t, x)  # 真值
        cen_difference()
    

    代码运行结果:

    16的二阶中心差分值:29.6942056593151
    16的二阶中心差分的误差:0.0691570640275347%
    

    二阶中心差分

    在这里插入图片描述
    Python代码:

    import sympy
    from sympy import diff
    from sympy import symbols
    
    #差分的对象
    x = 16
    k = 2  #步长
    x1 = x + k  #向前
    x2 = x - k  #向后
    
    #方程式
    def func(t):
        return 2000 * sympy.log(14*10000/(14*10000-2100*t))-9.8*t
    
    #二阶中心差分
    def two_cen_difference():
        a_cen_diff = (func(x1)+func(x2)-2*func(x))/(k**2)
        cen_error = abs(a_cen_diff - a_true)/a_true
        print(f'{x}的二阶中心差分值:{a_cen_diff}')
        print(f'{x}的二阶中心差分的误差:{cen_error*100}%')
    
    if __name__ == '__main__':
        t = symbols("t")
        a_true = diff(func(t), t , 2).subs(t, x)  # 求2阶导真值
        two_cen_difference()
    

    程序运行结果:

    16的二阶中心差分值:0.779693478624694
    16的二阶中心差分的误差:0.0779896119162236%
    

    迭代直到精确要求——以一阶向前差分为例

    要求初始步长为2,经过多次迭代,一阶向前差分的误差能小于1%
    Python代码如下:

    import sympy
    from sympy import diff
    from sympy import symbols
    
    #方程式
    def func(t):
        return 2000 * sympy.log(14*10000/(14*10000-2100*t))-9.8*t
    
    #一阶向前差分
    def for_difference(x1,x,a_true,k):
        a_for_diff = (func(x1) - func(x))/k
        for_error = abs(a_for_diff - a_true)/a_true
        print(f'{x}的一阶向前差分值:{a_for_diff}')
        print(f'{x}的一阶向前差分的误差:{for_error*100}%')
        return for_error
    
    def Judge_precision(x):
        D = True
        k = 2  # 初始步长
        n = 0
        while D:
            n += 1
            x1 = x + k  # 向前
            t = symbols("t")
            a_true = diff(func(t), t).subs(t, x)
            error = for_difference(x1,x,a_true,k)
            if error <= 0.01:
                D = False
                print(f'迭代第{n}次后,{x}的一阶向前差分的误差为{error}')
            else:
                k = k/2
    
    if __name__ == '__main__':
        x = 16  #差分对象
        Judge_precision(x)
    

    运行结果:

    16的一阶向前差分值:30.4738991379398
    16的一阶向前差分的误差:2.69671578943888%
    16的一阶向前差分值:30.0684298016344
    16的一阶向前差分的误差:1.33028844112341%
    16的一阶向前差分值:29.8697466293837
    16的一阶向前差分的误差:0.660728265039121%
    迭代第3次后,16的一阶向前差分的误差为0.00660728265039121
    
    展开全文
  • 二维Poisson方程
  • 偏微分方程的有限差分法的 MATLAB 和 Python 实现 这些代码实现了有限差分法的数值方法来求解 Heat PDE 和 Black-Scholes PDE。 具体而言,Black-Scholes PDE 的代码旨在为普通期权定价,例如欧洲和美国的看涨和看跌...

空空如也

空空如也

1 2 3 4 5
收藏数 95
精华内容 38
关键字:

python求解差分方程

python 订阅