Решение системы УЧП второго порядка с использованием Рунге Кутты на C

У меня возникла проблема с решением системы дифференциальных уравнений с использованием алгоритма Рунге-Кутта. До сих пор я переписал УЧП второго порядка в систему двух связанных уравнений, где

    f(L1,L2) = L2
    g(L1,L2) = A*(B*L1-C*L2-D)

- два уравнения, а A, B, C и D - константы. Чтобы получить значение для следующего шага, я проделал следующее для каждого временного шага dt:

    k1 = f(L1,L2)
    l1 = g(L1,L2)

    k2 = f(L1 + 0.5 * dt * k1,L2 + 0.5 * dt * l1 )
    l2 = g(L1 + 0.5 * dt * k1,L2 + 0.5 * dt * l1 )

    k3 = f(L1 + 0.5 * dt * k2,L2 + 0.5 * dt * l2 )
    l3 = g(L1 + 0.5 * dt * k2, L2 + 0.5 * dt * l2 )

    k4 = f(L1 + dt * k1,L2 +  dt * l1 )
    l4 = g(L1 + dt * k1,L2 + dt * l1 )

Где я использую значения для L1 и L2 текущего временного шага и итеративно вычисляю коэффициенты.

В результате я получаю L1 и L2, суммируя и взвешивая коэффициенты в конце. Моя проблема в том, что весь алгоритм становится нестабильным после 4 временных шагов.

Кто-нибудь знает, правильна ли реализация технически? Спасибо!


c system pde runge-kutta
person MichaelScott    schedule 08.11.2013    source источник
comment
Подскажите, пожалуйста, L1 и L2?   -  person Jonas Bötel    schedule 10.11.2013
comment
L1 и L2 - значения. для первого шага я использовал здесь начальные значения.   -  person MichaelScott    schedule 10.11.2013
comment
Под значениями вы имеете в виду константы, верно? Какие начальные значения вы используете? Можете ли вы показать результат вашей интеграции для пары шагов?   -  person Jonas Bötel    schedule 10.11.2013

Ответы (3)


arrow_upward
0
arrow_downward

просто предположение, так как вы не говорите, какое dt значение вы используете: сделайте его как можно меньше, потому что

локальная ошибка усечения имеет порядок O (h ^ 5), а общая накопленная ошибка порядка O (h ^ 4).

(цитируется из этой статьи в Википедии, dt играет роль h).

person CapelliC    schedule 10.11.2013
comment
Спасибо за Ваш ответ. Дело в том, что в моем случае это не зависит от dt (time inkrement). Если я положу dt = 0,001 или 0,01, большой разницы в значениях нет. Вот почему я думаю, что это неправильно ... У меня вопрос, правильно ли это реализовать Рунге Кутта, поскольку я понятия не имею, как это проверить по-другому ... - person MichaelScott; 10.11.2013

arrow_upward
0
arrow_downward

Две вещи:

Рунге-Кутта в целом нестабилен. Он просто «более устойчив», чем у Эйлера. В зависимости от условия вашего дифференциального уравнения и dt этого может быть просто недостаточно. Помогает ли меньший dt?

Мне не хватает понятия t в вашем определении f и g. Предполагая, что L1 и L2 не являются постоянными в t, вам лучше передать их через f и g. Вроде f(t,L1,L2). Это заставляет вас задуматься о тех расчетах коэффициентов, где вам теперь нужно передать еще t' соответственно. Это приведет к оценке L1 и L2 в средних точках.

person Jonas Bötel    schedule 10.11.2013

arrow_upward
0
arrow_downward

При вычислении k4, l4 вы использовали k1, l1 для вычисления состояния, в котором оценивается функция ODE, в то время как метод требует использования k3, l3. Эта ошибка уменьшит порядок метода, вероятно, до 2, и зависимость ошибки от размера шага будет отражать это.

person Lutz Lehmann    schedule 15.06.2018