Методы Рунге – Кутты для интеграции ОДУ в Python с дополнительными ограничениями - PullRequest
Новая
       135

Методы Рунге – Кутты для интеграции ОДУ в Python с дополнительными ограничениями

1 голос
/ 21 июня 2020

У меня вопрос о решении дифференциальных уравнений второго порядка с помощью RK4 с учетом дополнительных ограничений на первую производную. Я делаю показанный пример здесь с некоторыми изменениями

θ ′ ′ (t) + b θ ′ (t) + c sin (θ (t )) = 0

Дополнительное ограничение:

при θ i θ (i + 1) <0, тогда <em>θ ′ (i + 1) = 0,9 θ ′ i ,

где i - шаги t, i + 1 - один шаг после i. В реальном мире он говорит, что когда направление смещения меняется на противоположное, его скорость уменьшается до 90%.

Векторно, если y ( t ) = ( θ ( t ), ω ( t )), тогда уравнение будет ẏ = f ( t , y ), где f ( t , y ) = ( y ₂ ( t ), - by ₂ ( t ) - cos ( y ₁ ( t ))).

Как мне изменить код в этой задаче, если я хочу добавить ограничения на ω или θ ′ (t) (что одно и то же)? Вот мой код, который не сработал. Дополнительное условие делает θ не непрерывным. Следующее «самодельное» решение не может правильно обновить θ ′. Я новичок в Python, и это мой первый пост на StackOverflow. Приветствуются любые указания. 

def rungekutta4(f, y0, t, args=()):
    n = len(t)
    y = np.zeros((n, len(y0)))
    y[0] = y0
    for i in range(n - 1):
        h = t[i+1] - t[i]
        if y[i][0]*y[i+1][0]<0:
            k1 = f(y[i], t[i], *args)
            k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
            k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
            k4 = f(y[i] + k3 * h, t[i] + h, *args)
            y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)*0.9
        else:
            k1 = f(y[i], t[i], *args)
            k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
            k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
            k4 = f(y[i] + k3 * h, t[i] + h, *args)
            y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)        
    return y

Ответы [ 2 ]

1 голос
/ 21 июня 2020

Если я вас полностью не понимаю, вам нужно простое различение регистра в f: Математически вы получите f₂ ( t , y ) = - на ₂ ( t ) - cos ( y ₁ ( t )) если θ i ² − 1 = y ₁ ( t ) ² − 1 <0 и 0,9 · (<em> y ₂ − 1) в противном случае. Это все еще только зависимость f от y, , которая может быть просто реализована программно.

Например, вы можете определить: 

def f(y):
    θ = y[0]
    ω = y[1]
    return [
        θ,
        -b*ω-cos(θ) if θ**2<1 else 0.9*(ω-1)
    ]

Хотя это может работать как есть, вы можете столкнуться с проблемами из-за того, что f не является непрерывным (или имеет непрерывную производную), в частности, если вы хотите использовать более продвинутые интеграторы с размером шага контроль вместо домашнего. Чтобы избежать этого, замените конструкцию if - else на (острый) сигмоид. Подробнее об этом см. этот мой ответ .

0 голосов
/ 22 июня 2020

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

    for i in range(n - 1):
        h = t[i+1] - t[i]
        y[i+1] = RK4step(f,t[i],y[i],h, args)
        if y[i+1,0]*y[i,0] < 0: y[i+1,1] *= 0.9
    return y

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

...