用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

使用： ^{cd5}

二阶导数方程的初始条件是p(0) = 10^14，p0 = 10^13，L = 10 ^-6和{}，第二个条件可能会出错。在

我试着用直接的方法解决这个问题-类似于前面的

^{pr2}$

但结果是不正确的。它应该给我一个指数排序的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

用法示例：

^{pr2}$

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