Дифференциальные уравнения

Рейтинг: 0Ответов: 2Опубликовано: 06.07.2023

Здравствуйте у меня есть функция ускорения следующего вида a(x, t), 2 начальный условия x(0) = x0, v(0) = v0. Нужно восстановить функции v(t), x(t) при 0< t <t0. Я пробовал следующее: разбивал время на малые сегменты dt =0.01, затем считал , что на малом сегменты движение - равноускоренное, т. е. Δx = v *dt + a * dt * dt / 2, Δv = a * dt. Такие рассуждения привели к неверным результатам. Также хотел использовать функции для решения дифференциальных уравнений, но т. к. a(x, t) - функция 2 переменных, то ничего не получалось.

Скажите пожалуйста, есть ли какие - нибудь функции в python для численного решения такого дифференциального уравнения?

Ответы

▲ 0Принят

В пакете scipy.integrate есть целая пачка функций для интегрирования обыкновенных дифференциальных уравнений: https://docs.scipy.org/doc/scipy/reference/integrate.html#solving-initial-value-problems-for-ode-systems

Универсальная обёртка - функция solve_ivp

Пример: предмет подбрасывают вверх, начальная высота 10 метров, скорость 1 метр в секунду. Решение вычисляется для интервала t \in [0,3], используется метод по-умолчанию явный Рунге-Кутты пятого порядка

import scipy.integrate as spi
import numpy as np


def a(x, t):
    # Притяжение не зависит ни от высоты, ни от времени
    return -9.81


def diff_equation(t, xv):
    # Первый параметр функции - время, второй - вектор [высота, скорость]
    x, v = xv
    # dx/dt = v, dv/dt = a(x,t)
    return [v, a(x, t)]


x0 = 10
v0 = 1
# Вызов решателя, интервал 0-3 сек, начальные условия, набор 10 точек t от 0 до 3
solution = spi.solve_ivp(diff_equation, [0, 3], [x0, v0], t_eval=np.linspace(0,3, 10))
if solution.success:
    print("t: ", solution.t)
    print("x: ", solution.y[0])
    print("v: ", solution.y[1])
else:
    print("НЕШМОГЛА... ", solution)

Результат

t:  [0.         0.33333333 0.66666667 1.         1.33333333 1.66666667
 2.         2.33333333 2.66666667 3.        ]
x:  [ 10.           9.78833333   8.48666667   6.095        2.61333333
  -1.95833333  -7.62       -14.37166667 -22.21333333 -31.145     ]
v:  [  1.    -2.27  -5.54  -8.81 -12.08 -15.35 -18.62 -21.89 -25.16 -28.43]

Если 10 отсечек вам мало, то поменяйте значение параметра t_eval Например, t_eval=np.arange(0,3.001, 0.01) задаёт набор моментов времени от 0 до 3 секунд с шагом 0.01 секунда.

Набор точек в t_eval не влияет на точность вычислений, это лишь команда решателю для каких отсечек вывести значения.

▲ 0

Предварительные замечания на ваше

разбивал время на малые сегменты dt =0.01, затем считал , что на малом сегменты движение - равноускоренное, т. е. Δx = v *dt + a * dt * dt / 2, Δv = a * dt.

Идея примерно та, что надо, член второго порядка можно выбросить, оставив

dx/dt = v
dv/dt = a

И тогда получается на каждом участке просчет как

v(t+dt) = v(t) + a(x,t)*dt
x(t+dt) = x(t) + v(t)*dt

Получается простейшая схема Эйлера. Дальше зависит от требуемой точности, вида функции etc etc. Может оказаться достаточно этого, может, потребуется более точная схема Рунге-Кутты, может, что еще...

Вы утверждаете, что приведенная схема у вас не работает - как именно она не работает? Может, вы ее неверно реализовали? Без конкретного кода, вида функции и т.д. говорить не о чем. При нечетком ТЗ результат один — ХЗ...