2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте его в существующую тему, а создайте новую в корневом разделе "Помогите решить/разобраться (М)".

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

Не ищите на этом форуме халяву, правила запрещают участникам публиковать готовые решения стандартных учебных задач. Автор вопроса обязан привести свои попытки решения и указать конкретные затруднения.

Обязательно просмотрите тему Правила данного раздела, иначе Ваша тема может быть удалена или перемещена в Карантин, а Вы так и не узнаете, почему.



Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение01.12.2013, 12:46 


15/05/10
20
Дана система диф. уравнений:

$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 порядка (решения данной системы лишь подзадача). Однако результаты неправильные.

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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение01.12.2013, 21:48 


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение01.12.2013, 23:33 
Заслуженный участник


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение05.12.2013, 19:34 


15/05/10
20
Vince Diesel в сообщении #795206 писал(а):
Есть книга Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. Там много чего должно быть. Насколько помню, также прилагаются листинги готовых програм на фортране.


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

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

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

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


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение06.12.2013, 23:44 


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

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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение07.12.2013, 20:09 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение08.12.2013, 12:07 


15/05/10
20
10110111 в сообщении #797161 писал(а):
Цитата:
Использовал разные шаги, при больших например = 1, или = 10, результат уходит в NaN, при маленьких ответ не верен

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


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

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

Утундрий в сообщении #797456 писал(а):
Обычно в таких делах подключают вручную вычисленную матрицу Якоби.


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение08.12.2013, 14:44 
Заслуженный участник
Аватара пользователя


15/10/08
10065
Antistas в сообщении #797644 писал(а):
вычислить ее вручную, а затем вставить в программу в виде отдельной функции?

Натюрлихъ.

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение08.12.2013, 17:52 


15/05/10
20
Утундрий в сообщении #797708 писал(а):
Antistas в сообщении #797644 писал(а):
вычислить ее вручную, а затем вставить в программу в виде отдельной функции?

Натюрлихъ.


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение21.06.2014, 17:12 
Заслуженный участник


15/05/05
3444
USA
Antistas в сообщении #794889 писал(а):
Какие методы помимо Рунге Кутта можно использовать для решения подобных систем ?
Метод Адамса.

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение08.08.2014, 07:35 


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение08.08.2014, 09:24 
Заслуженный участник
Аватара пользователя


01/08/06
2868
Уфа
Метод Гира уже советовали?

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение21.01.2016, 17:43 


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение13.10.2018, 15:23 
Аватара пользователя


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

 Профиль  
                  
 
 Re: Неявный метод Рунге Кутта для решения жесткой системы ДУ
Сообщение10.08.2020, 15:14 
Аватара пользователя


30/04/19
209
Тема не новая, но вроде как в верхней части списка, поэтому отвечу. Как-то писал самый простой неявный метод Эйлера, правда для решения одного жесткого ДУ. Очень долго я его не тестировал, но на жестких уравнениях проверил, это хорошо работало там где явные схемы шли в разнос. В качестве решателя нелинейного уравнения использован метод Ньютона, в качестве начального приближения используется 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 Почему то отступы в коде удваиваются, становятся шире чем в среде разработки.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 18 ]  На страницу 1, 2  След.

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group