python - 使用 python 进行二阶数值微分

Question

我想使用Python计算二阶微分方程而不使用内置函数，但结果仅对于一阶方程是正确的。

让我给你举个例子(插图没有声誉! - 为了更好的方程外观)

dy/dt = -ky

并使用导数的基本定义

f'(x)= h ->0 (f(x+h)-f(x))/h

我们可以为 k=0.3 的方程编写基本的 Python 代码

def first_order(dt):
    t = np.arange(0, 20, dt)
    k = 0.3
    y = np.zeros(len(t))

    y[0] = 5
    for i in range(1, len(t)):
        y[i] = - k* y[i - 1] * dt + y[i-1]

    return t, y

这工作得很好，但是当我尝试计算相应的方程时:

dp^2/dx^2 = (p- p0)/L

使用: f''(x)= h ->0 (f(x+h)-2f(x)+ f(x-h)/h^2

二阶导数方程的初始条件为p(0) = 10^14、p0 = 10^13、L = 10 ^-6 和 p(infinity) = p0，第二个条件可能会导致错误。

我尝试以简单的方式解决这个问题 - 类似于之前的

def diffusion_lenght(dt):
   p0 = 10 ** 13  # initial state
   t = np.arange(0, 20, dt)
   p = np.zeros(len(t))
   L = 1 * 10 ** -6
   p[0] = 10 ** 14
   p[len(t)-1] = p0

   for i in range(1, len(t)):
      p[i] = (2* p[i-1]- p[i-2]- p0 * dt ** 2 / L) / (1 - dt ** 2 / L)
   print(p)
   return t, p

但结果不正确。它应该给我 x 的指数递减，但我得到了收敛到 dt 值的直线。

Answer 1

如果您不介意，我几年前就这样做了，并且我有一个小脚本，所以我会给您我的代码，而不是在您的代码中查找错误。

我认为它非常清晰明确。

import numpy as np

class dif_eq(object):
    def __init__(self):
        pass
    def ft4(self,dt,u,x1_before,x2_before,x3_before,x4_before,functionn):
        x1_now = x1_before + dt * x2_before
        x2_now = x2_before + dt * x3_before
        x3_now = x3_before + dt * x4_before
        x4_now = x4_before + dt * functionn(u,x1_before,x2_before,x3_before,x4_before)
        vals = [x1_now,x2_now,x3_now,x4_now]
        return vals

    def ft3(self,dt,u,x1_before,x2_before,x3_before,functionn):
        x1_now = x1_before + dt * x2_before
        x2_now = x2_before + dt * x3_before
        x3_now = x3_before + dt*functionn(u,x1_before,x2_before,x3_before)
        vals = [x1_now,x2_now,x3_now]
        return vals

    def ft2(self,dt,u,x1_before,x2_before,functionn):
        x1_now = x1_before + dt * x2_before
        x2_now = x2_before+  dt * functionn(u,x1_before,x2_before)
        vals = [x1_now,x2_now]
        return vals

    def ft1(self,dt,u,x1_before,functionn):
        x1_now = x1_before + dt*functionn(u,x1_before)
        vals = [x1_now]
        return vals

使用示例:

"""
def order3equation(u,b,c,a): # Your differential equations
    y=u-b*2-c*2.5-a*3.6 #+ noise*np.random.rand()/5
    return y
def order1equation.....
    ....
    return y

d=dif_eq()
val=[0] # init for 1order
val3=[0,0,0] # init for 3order
dt=0.05
result1order=[]
result3order=[]
for i in range(100):
    val=d.ft1(dt,u,val[0],order1equation)
    result1order.append(val[0])

for i in range(1000):
    val3=d.ft3(dt,u,val3[0],val3[1],val3[2],order3equation)
    result3order.append(val3[0])
"""

valx/vals 是导数的实际值。第一次传递初始值，然后将函数实际返回的值传递给函数。

工作示例 - 不稳定的系统

    u = 1 # input/impulse or whatever name it is
    d=dif_eq()

    val3=[0,0,0]
    dt=0.05

    zz=[]
    def my3(u, b, c, a):
        y = u - b * 1 - c * 1 - a *1
        return y
    for i in range(1000):
        val3=d.ft3(dt,u,val3[0],val3[1],val3[2],my3)
        zz.append(val3[0])
    from matplotlib.pyplot import *
    plot(zz)
    show()