Неявный метод Рунге Кутта для решения жесткой системы ДУ

[Antistas]

Дана система диф. уравнений:

$c_1 ' = (-k_3 c_1 (1 - x_1 - x_2 - x_3 - x_4) + k_3 x_3) \cdot L / N + (c_1 ^ {0} - c_1) \cdot V / N$ $c_2 ' = (-k_1 c_2 (1 - x_1 - x_2 - x_3 - x_4) + k_{-1} x_1 - k_2 c_2 (1 - x_1 - x_2 - x_3 - x_4) ^ {2} + k_{-2} x_2 ^ 2) \cdot L / N + (c_2 ^ {0} - c_2) \cdot V / N$ $c_3 ' = k_5 x_1 x_4 L / N + (c_3^0 - c_3) V / N$ $c_4 ' = k_5 x_1 x_4 L / N + (c_4 ^ 0 - c_4) V / N$ $x_1 ' = k_1 c_2 (1 - x_1 - x_2 - x_3 - x_4) - k_{-1} x_1 - k_5 x_1 x_4$ $x_2 ' = 2(k_2 (1 - x_1 - x_2 - x_3 - x_4) ^ 2 c_2 - k_{-2} x_2 ^ 2) - k_4 x_2 x_3$ $x_3 ' = k_3(1- x_1 - x_2 - x_3 - x_4) c_1 - k_3 x_3 - k_4 x_2 x_3$ $x_4 ' = 2 k_4 x_2 x_3 - k_5 x_1 x_4$ $T' = (T^0 - T) V / N + Q k_5 x_1 x_4 / C_P$

Константы $k_{-3}, k_{-2}, k_{-1}, k_{1}, k_{2}, k_{3}, k_{4}, k_{5}$ разных порядков, примерно от $10^{11}$ до $10^{24}$. Поэтому система является довольно жесткой. Пытаюсь решить неявным метом - с целью простоты методом Рунге Кутта 2 порядка (решения данной системы лишь подзадача). Однако результаты неправильные.

Какие методы помимо Рунге Кутта можно использовать для решения подобных систем ?

[sithif]

На сколько я понимаю, неявный метод не гарантирует "правильного" решения.
Для неявных схем просто возможно использовать больший шаг по сравнению с явными.
Если у вас нет ошибок в самом алгоритме, попробуйте уменьшать шаг интегрирования.
Можете попробовать найти ограничение на шаг интегрирования спектральным методом.
В мануале Mathematica есть, что то про жесткие системы, также там ссылки на какие-то статьи по теме.
http://www.wolfram.com/learningcenter/tutorialcollection/AdvancedNumericalDifferentialEquationSolvingInMathematica/

[Vince Diesel]

Есть книга Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. Там много чего должно быть. Насколько помню, также прилагаются листинги готовых програм на фортране.

[Antistas]

Vince Diesel в сообщении #795206 писал(а):

Есть книга Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. Там много чего должно быть. Насколько помню, также прилагаются листинги готовых програм на фортране.

Да, спасибо, я это все видел. Я пишу программу сам, а вот метод который все это дело решает не получается.. Или такая вот система.

Могут ли здесь на форуме помочь с программой, если я скину исходный код ?

-- Чт дек 05, 2013 20:35:51 --

sithif в сообщении #795185 писал(а):

На сколько я понимаю, неявный метод не гарантирует "правильного" решения.
Для неявных схем просто возможно использовать больший шаг по сравнению с явными.
Если у вас нет ошибок в самом алгоритме, попробуйте уменьшать шаг интегрирования.
Можете попробовать найти ограничение на шаг интегрирования спектральным методом.
В мануале Mathematica есть, что то про жесткие системы, также там ссылки на какие-то статьи по теме.
http://www.wolfram.com/learningcenter/tutorialcollection/AdvancedNumericalDifferentialEquationSolvingInMathematica/

Использовал разные шаги, при больших например = 1, или = 10, результат уходит в NaN, при маленьких ответ не верен

[10110111]

Цитата:

Использовал разные шаги, при больших например = 1, или = 10, результат уходит в NaN, при маленьких ответ не верен

Вероятно, сказывается ошибка округления. Попробуйте использовать более длинные представления чисел (например, long double вместо double, или вообще длинную арифметику).

[Утундрий]

Обычно в таких делах подключают вручную вычисленную матрицу Якоби.

[Antistas]

10110111 в сообщении #797161 писал(а):

Цитата:

Использовал разные шаги, при больших например = 1, или = 10, результат уходит в NaN, при маленьких ответ не верен

Вероятно, сказывается ошибка округления. Попробуйте использовать более длинные представления чисел (например, long double вместо double, или вообще длинную арифметику).

Использую Extended вроде как этого хватает, числа переваливают за E+400

-- Вс дек 08, 2013 13:08:52 --

Утундрий в сообщении #797456 писал(а):

Обычно в таких делах подключают вручную вычисленную матрицу Якоби.

Ну все таки мне ее придется считать на каждом шаге , т.е. раз 200. Или вы мне предлагаете вычислить ее вручную, а затем вставить в программу в виде отдельной функции?

[Утундрий]

Antistas в сообщении #797644 писал(а):

вычислить ее вручную, а затем вставить в программу в виде отдельной функции?

Натюрлихъ.

[Antistas]

Утундрий в сообщении #797708 писал(а):

Antistas в сообщении #797644 писал(а):

вычислить ее вручную, а затем вставить в программу в виде отдельной функции?

Натюрлихъ.

Легко сказать:) для метода второго порядка придется только 18 * 18 функций писать:)

[Yuri Gendelman]

Antistas в сообщении #794889 писал(а):

Какие методы помимо Рунге Кутта можно использовать для решения подобных систем ?

Метод Адамса.

[Glk63]

Я использую метод Розенброка с действительными или комплексными коэффициентами. С комплексными лучше, но считает существенно дольше.

[worm2]

Метод Гира уже советовали?

[saygogoplz]

А напишите константы в численном виде, можно попробывать сравнить результаты (мои и ваши).

[Eugeen1948]

Согласно теории "жесткость" ОДУ определяется соотношением наибольшего собственного числа (СЧ) матрицы Якоби правых частей ОДУ к наименьшему. Величина констант, входящих в уравнения может только "намекать" на жесткость.
Еще сложнее случай, когда спектр СЧ комплексный. Многие методы в этом случае "буксуют".
Есть "подводные камни" при автоматическом выборе шага интегрирования (или порядка точности метода, напр. у метода Гира). Это характерно для решений, которые либо пересекают 0 либо колеблются около 0. Тогда для контроля шага/порядка надо искусственно вводить масштаб погрешности для этих переменных.
В настоящее время пользуются разными пакетами-решателями, основанными на методе BDF (см. в инете). Собственно метод Гира и есть BDF с фиксированным шагом, где пересчет шага можно делать только по прошествии нескольких шагов с постоянным шагом. Для "чистого" BDF это ограничение снято.
Мне представляется, что для экстремально жестких систем ДУ надо применить более центрированную формулу "дифференцирования назад" добавив кроме неявного слоя (n+1) слой (n+2). Это один из вариантов т.н. "Неявных Блочных Методов", предложенный Клиппинджером и Димсдейлом (см. Холл Дж. Уатт Дж. (ред.) Современные численные методы решения обыкновенных дифференциальных уравнений, Мир, 1979. - 312с.).
Вообще лучше книги по численным методам для ОДУ я не встречал, хотя в моей библиотеке таких более 100. Вы точно найдете в ней полезные рекомендации по написанию своей программы.

[Snegovik]

Тема не новая, но вроде как в верхней части списка, поэтому отвечу. Как-то писал самый простой неявный метод Эйлера, правда для решения одного жесткого ДУ. Очень долго я его не тестировал, но на жестких уравнениях проверил, это хорошо работало там где явные схемы шли в разнос. В качестве решателя нелинейного уравнения использован метод Ньютона, в качестве начального приближения используется Y найденный на предыдущем шаге интегрирования, никаких дополнительных "сложностей" я не вводил. Единтвенное до чего не добрались руки - протестить вдоль и поперек и установить правильные значения eps (в методе Ньютона) и dy (там же, для численного поиска производной).

Используется синтаксис Fortran

!Модуль решения диф уравнения неявным методом Эйлера
module eyler
implicit none

public :: eyler_n
private 

        contains

!Подпрограмма решения диф. уравнения
subroutine eyler_n(x0, x1, n, y0, y1, error)
real (8), intent (in)   :: x0           !Начальное Х
real (8), intent (in)   :: x1           !Конечное Х
real (8), intent (in)   :: y0           !Начальное Y
real (8), intent (out)  :: y1           !Искомое Y
real (8)                :: x            !Текущее X
real (8)                :: y            !Текущее Y       
real (8)                :: h            !Шаг интегрирования
integer                 :: n            !Число шагов интегрирования
integer                 :: i            !Счетчик цикла
integer, intent (out)   :: error        !Сообщение об ошибке

if (n < 30) n = 30
h = (x1-x0) / float(n)
x = x0 + h
y = y0

        do i = 1, n
                call newton(x, y, h, y, y1, error)  
                y = y1
                x = x + h 
                if (error == 1) exit
        end do

end subroutine eyler_n

!Подпрограмма решения нелинейного уравнения методом Ньютона
subroutine newton (x, y, h, y0, y2, error)
real (8), intent (in)   :: x                            !Текущее Хi+1
real (8), intent (in)   :: y                            !Текущее Yi
real (8), intent (in)   :: h                            !Шаг интегрирования
real (8), intent (in)   :: y0                           !Начальное приближение для Yi+1
real (8), intent (out)  :: y2                           !Искомое значение Yi+1
real (8), parameter     :: dy = 1.0e-6                  !Дифференциал по Y
real (8), parameter     :: eps = 1.0e-5                 !Необходимая точность
real (8)                :: ytemp1, ytemp2               !Временные Y
real (8)                :: dg                           !Производная функции g        
integer                 :: i, error                     !Счетчик и сообщение о зацикливании

        error = 0
        ytemp1 = y0
        i=0
        
                do while (abs(g(y, h, x, ytemp2)) > eps)
                        dg = (g(y, h, x, ytemp1+dy) - g(y, h, x, ytemp1-dy)) / (dy+dy)
                        ytemp2 = y - g(y, h, x, ytemp1) / dg
                        ytemp1 = ytemp2
                        i = i + 1
                          if (i > 10000) then
                                error = 1
                                exit
                          end if
                end do
        
        y2 = ytemp2             

end subroutine newton
        
 !Пробразованная функция для решения методом Ньютона
 pure real (8) function g(y, h, x, y1)
 real (8), intent(in) :: y      !Текущее Y
 real (8), intent(in) :: h      !Шаг
 real (8), intent(in) :: x      !Текущее Х
 real (8), intent(in) :: y1     !Подставляемое Y1
                
        !Интерфейс к функции дифференциального уравнения
        interface
          pure real (8) function f (x, y)
            real (8), intent (in) :: x, y
          end function f
        end interface   
        
 g = y + h * f(x, y1) - y1
 end function g
        
end module eyler

 

Тестирующая программа и функция диф. уравнения могут выглядеть так.

Используется синтаксис Fortran

!Главная программа
program test
use eyler
implicit none
real (8)        :: x0           !Начальное X 
real (8)        :: x1           !Конечное X
real (8)        :: y0           !Начальное Y
real (8)        :: y1           !Конечное Y
integer         :: n            !Число шагов интегрирования
integer         :: error        !Обшибка метода Ньютона (зацикливание)

!Начальные условия
x0  = 0.0
x1  = 3.0
y0  = 1.0     
n   = 500

 call eyler_n (x0, x1, n, y0, y1, error)

print '(a, 1p, e15.7)', "Решение   = ",  y1
print '(a, i3)', "Ошибка           = ",  error
end program test

!Функция диф. уравнения
pure real (8) function f (x, y)
        real (8), intent (in) :: x, y
        f = 50 * (exp(-x)-y)
end function f

 

Заодно интересно, вдруг кто обнаружит ошибку (если есть конечно).

ps Почему то отступы в коде удваиваются, становятся шире чем в среде разработки.