Дифференциал второго порядка численно с питоном - PullRequest
Новая
       28

Дифференциал второго порядка численно с питоном

0 голосов
/ 19 октября 2018

Я хотел бы рассчитать дифференциальное уравнение второго порядка, используя Python без использования встроенных функций, но результаты верны только для уравнения первого порядка.

Позвольте мне привести вам пример (Нет репутации для врезанных изображений! - для лучшего вида уравнения)

dy/dt = -ky

И используя базовое определение производной

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

Мы можем написать основной код Python для этого уравнения для k = 0,3 

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.

1 Ответ

0 голосов
/ 19 октября 2018

Если вы не возражаете, я сделал это несколько лет назад, и у меня есть небольшой скрипт для него, поэтому я дам вам свой код вместо того, чтобы находить ошибки в вашем.

Я думаю, что он очень ясный и явный.

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