![]() |
|
сделать стартовой | добавить в избранное |
![]() |
Решение систем дифференциальных уравнений методом Рунге-Куты 4 порядка |
Министерство образования Украины Донецкий государственный технический университет Кафедра химической технологии топлива Курсовая работана тему : Решение систем дифференциальных уравнений методом Рунге - Кутты 4 порядкапо дисциплине : Математические методы и модели в расчетах на ЭВМВыполнил: студент гр. ХТ-96 Кузнецов М.В.Проверил: доц. Чеховской Б.Я. г. Донецк 1998 год РЕФЕРАТ Дифференциальные Уравнения, Метод Рунге-Кутта, РК-4, Концентрация, Метод Эйлера, Задача Коши, Ряд Тейлора, Паскаль, Реакция, Интервал, Коэффициенты Дифференциального Уравнения.Листов : 28 Таблиц : 2 Графиков : 4 Решить систему дифференциальных уравнений методом Рунге-Кутты 4 порядка, расчитать записимость концентрации веществ в зависимости от времени, проанализировать полученную зависимость, удостовериться в действенности метода. Содержание: Введение1. Постановка задачи 62. Суть метода 83. Выбор метода реализации программы 144. Блок – схема .155. Программа .176. Идентификация переменных 197. Результаты .208. Обсуждение результатов .219. Инструкция к программе .2310. Заключение .27 Литература Введение Обыкновенные дифференциальные уравнения (ОДУ) широко используются для математического моделирования процессов и явлений в различных областях науки и техники. Переходные процессы в радиотехнике, кинетика химических реакций, динамика биологических популяций, движение космических объектов, модели экономического развития исследуются с помощью ОДУ. В дифференциальное уравнение -го порядка в качестве неизвестных величин входят функция y(x) и ее первые производных по аргументу x (( x, y, y1, . y( ) )=0. 1.1 Из теории ОДУ известно, что уравнение (1.1) эквивалентно системе уравнений первого порядка (k(x, y1, y1’ ,y2 ,y2 ’, . ,y ,y ’)=0. 1.2где k=1, . , . Уравнение (1.1) и эквивалентная ему система (1.2) имеют бесконечное множество решений. Единственные решения выделяют с помощью дополнительных условий, которым должны удовлетворять искомые решения. В зависимости от вида таких условий рассматривают три типа задач, для которых доказано существование и единственность решений. Первый тип – это задачи Коши, или задачи с начальными условиями. Для таких задач кроме исходного уравнения (1.1) в некоторой точке xo должны быть заданы начальные условия, т.е. значения функции y(x) и ее производных y(x0)=y0’ , y’(x0)=y10, . , y( -1)(x0)=y -1,0. Для системы ОДУ типа (1.2) начальные условия задаются в виде y1(x0)=y10 , y2(x0)=y20, . , y (x0)=y 0. 1.3 Ко второму типу задач относятся так называемые граничные, или краевые задачи, в которых дополнительные условия задаются в виде функциональных соотношений между искомыми решениями. Количество условий должно совпадать с порядком уравнения или системы. Если решение задачи определяется в интервале x є (x0 ,xk(, то такие условия могут быть заданы как на границах, так и внутри интервала. Минимальный порядок ОДУ, для которых может быть сформулирована граничная задача, равен двум. Третий тип задач для ОДУ – это задачи на собственные значения. Такие задачи отличаются тем, что кроме искомых функций y(x) и их производных в уравнения входят дополнительно m неизвестных параметров (1((2((( хm( которые называются собственными значениями( Для единственности решения на интервале необходимо задать m граничных условий( В качестве примера можно назвать задачи определения собственных частот( коэффициентов диссипации( структуры электромагнитных полей и механических напряжений в колебательных системах( задачи нахождения фазовых коэффициентов( коэффициентов затухания( распределения напряженностей полей волновых процессов и т(д( К численному решению ОДУ приходится обращаться( когда не удается построить аналитическое решение задачи через известные функции( Хотя для некоторых задач численные методы оказываются более эффективными даже при наличии аналитических решений( Большинство методов решения ОДУ основано на задаче Коши( алгоритмы и программы для которой рассматриваются в дальнейшем( 1.
Постановка задачи Многие процессы химической технологии описываются СДУ - начиная от кинетических исследований и заканчивая химическими технологическими процессами( В основу математических способов описания процессов положены СДУ и СЛАУ( Эти уравнения описывают материальные и тепловые балансы объектов химической технологии( а так же структуры потоков технических веществ в этих аппаратах( Для получения( распределения технологических параметров во времени и в пространстве (в пределах объекта)( необходимо произвести СДУ методом( которых дал бы высокую точность решения при минималььных затратах времени на решение( потому что ЭВМ должна работать в режиме реального времени и успевать за ходом технологического процесса( Если время на решение задачи большое( то управляющее воздействие( выработанное на ЭВМ может привести к отрицательным воздействиям( Методов решения существует очень много( В данной работе будет рассмотрен метод решения СДУ методом Рунге-Кутта 4 порядка. Для удобства работы на ЭВМ, необходимо данную кинетическую схему преобразовать в удобный для работы на компьютере вид. Для этого необходимо кинетическую схему процесса представить в виде уравнений. При рассмотрении кинетической схемы процесса необходимо учитывать коэффициенты скоростей реакций. Но, так как процесс протекает при изотермических условиях, коэффициенты скоростей реакций можно считать за константы скоростей химической реакции. Из приведенной ниже схемы мы можем составить ряд дифференциальных уравнений, учитывающих изотермичность процесса.Так как коэффициенты K1,K2,K3,K4 являются константами, то можно уравнение записать в следущем виде. Для преобразования данных дифференциальных уравнений для использования их в расчетах тепловых и кинетических схем методами Рунге-Кутты необходимо подставлять вместо производных значений концентраций, значения концентраций данных в начале процесса. Это обусловлено тем, что метод Рунге-Кутты четвертого порядка, который будет использован для расчета кинетической схемы процесса. Так как этот метод требует сведений только об одной точке и значений функции. 2. Суть метода Разбор и рассмотрение методов( применяемых на практике для решения дифференциальных уравнений( мы начнем с их широкой категории( известной под общим названием методов Рунге-Кутта( Методы Рунге-Кутта обладают следующими свойствами:1( Эти методы являются одноступенчатыми: чтобы найти уm 1( нужна информация о предыдущей точке xm(ym(2( Они согласуются с рядом Тейлора вплоть до членов порядка hp( где степень р различна для различных методов и называется порядковым номером или порядком метода(3( Они не требуют вычисления производных от f (x(y)( а требуют вычисления самой функции( Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий( После этого мы подтвердим полученные результаты аналитически( Предположим( нам известна точка xm(ym на искомой кривой( Тогда мы можем провести прямую линию с тангенсом угла наклона у(m=f(xm(ym)( которая пройдет через точку xm(ym( Это построение показано на рис(1( где кривая представляет собой точное( но конечно неизвестное решение уравнения( а прямая линия L1 построена так( как это только что описано(Тогда следующей точкой решения можно считать ту( где прямая L1 пересечет ординату( проведенную через точку x=xm 1=xm h(Уравнение прямой L1 выглядит так: y=ym y(m(x-xm) так как y(=f(xm(ym) и кроме того( xm 1=xm h тогда уравнение примет вид ym 1=ym h f(xm(ym) 1(1 Ошибка при x=xm 1 показана в виде отрезка е( Очевидно( найденное таким образом приближенное значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h( так что ошибка ограничения равна e =Кh2 Заметим( что хотя точка на графике 1 была показана на кривой( в действительности ym является приближенным значением и не лежит точно на кривой( Формула 1(1 описывает метод Эйлера( один из самых старых и широко известных методов численного интегрирования дифференциальных уравнений( Отметим( что метод Эйлера является одним из методов Рунге-Кутта первого порядка( Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера( В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: xm(ym и xm h(ym hy(m( Последняя точка есть та самая( которая в методе Эйлера обозначалась xm 1(ym 1( Геометрический процесс нахождения точки xm 1(ym 1 можно проследить по рис(2( С помощью метода Эйлера находится точка xm h(ym hy(m( лежащая на прямой L1( В этой точке снова вычисляется тангенс( дает прямую (( Наконец( через точку xm(ym мы проводим прямую L( параллельную (( Точка( в которой прямая L пересечется с ординатой( восстановленной из x=xm 1=xm h( и будет искомой точкой xm 1(ym 1(Тангенс угла наклона прямой ( и прямой L равен Ф(xm(ym(h)=( 1(2 где y(m=f(xm(ym) 1(3 Уравнение линии L при этом записывается в виде y=ym (x-xm)Ф(xm(ym(h)( так что ym 1=ym hФ(xm(ym(h)( 1(4 Соотношения 1(2( 1(3( 1(4 описывают исправленный метод Эйлера( Чтобы выяснить( насколько хорошо этот метод согласуется с разложением в ряд Тейлора( вспомним( что разложение в ряд функции f(x(y) можно записать следующим образом: f(x(y)=f(xm(ym) (x-xm)(f/(x (y-ym)(f/(x ( 1(5где частные производные вычисляются при x=xm и y=ym( Подставляя в формулу 1(5 x=xm h и y=ym hy(m и используя выражение 1(3 для y(m( получаем f(xm h(ym hy(m)=f hfx hffy O(h2)(где снова функция f и ее производные вычисляются в точке xm(ym( Подставляя результат в 1(2 и производя необходимые преобразования( получаем Ф(xm(ym(h)=f h/2(fx ffy) O(h2)( Подставим полученное выражение в 1(4 и сравним с рядом Тейлора ym 1=ym hf h2/2(fx ffy) O(h3)( Как видим( исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h2( являясь( таким образом( методом Рунге- Кутты второго порядка( Рассмотрим модификационный метод Эйлера( Рассмотрим рис(3 где первоначальное построение сделано так же( как и на рис(2( Но на этот раз мы берем точку( лежащую на пересечении этой прямой и ординатой x=x h/2( На рисунке эта точка образована через Р( а ее ордината равна y=ym (h/2)y(m( Вычислим тангенс угла наклона касательной в этой точке Ф(xm(ym(h)=f (xm h/2(ym h/2 y(m)( 1(6 где y(m=f(xm(ym) 1(7 Прямая с таким наклоном( проходящая через Р( обозначена через ( ( Вслед за тем( мы проводим через точку xm(ym прямую параллельную ( ( и обозначаем ее через L0( Пересечение этой прямой с ординатой x=xm h и даст искомую точку xm 1(ym 1( Уравнение прямой можно записать в виде y=ym (x- xm)Ф(xm(ym(h)( где Ф задается формулой 1(6( Поэтому ym 1=ym hФ(xm(ym(h) 1(8 Соотношения 1(6( 1(7( 1(8 описывают так называемый модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго порядка( Обобщим оба метода( Заметим( что оба метода описываются формулами вида ym 1=ym hФ(xm(ym(h) 1(9и в обоих случаях Ф имеет вид Ф(xm(ym(h)=a1f(xm(ym) a2f(xm b1h(ym b2hy(m)( 1(10 где y(m=f(xm(ym) 1(11 В частности( для исправленного метода Эйлера a1=a2=1/2; b1=b2=1(В то время как для модификационного метода Эйлера a1=0( a2=1( b1=b2=1/2(Формулы 1(9( 1(10( 1(11 описывают некоторый метод типа Рунге-Кутты( Посмотрим( какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a1( a2( b1 и b2 ( Чтобы получить соответствие ряду Тейлора вплоть до членов степени h( в общем случае достаточно одного параметра( Чтобы получить согласование вплоть до членов степени h2( потребуется еще два параметра( так как необходимо учитывать члены h2fx и h2ffy( Так как у нас имеется всего четыре параметра( три из которых потребуются для создания согласования с рядом Тейлора вплоть до членов порядка h2( то самое лучшее( на что здесь можно рассчитывать - это метод второго порядка(В разложении f(x(y) в ряд 1(5 в окрестности точки xm(ym положим x=xm b1h( y=ym b2hf( Тогда f(xm b1h(ym b2hf)=f b1hfx b2hffy O(h2)( где функция и производные в правой части равенства вычислены в точке xm(ym(Тогда 1(9 можно переписать в виде ym 1=ym h O(h3)( Сравнив эту формулу с разложением в ряд Тейлора( можно переписать в виде ym 1=ym h O(h3)(Если потребовать совпадения членов hf( то a1 a2=1( Сравнивая члены( содержащие h2fx( получаем a2b1=1/2( Сравнивая члены( содержащие h2ffy( получаем a2b2=1/2(Так как мы пришли к трем уравнениям для определения четырех неизвестных( то одно из этих неизвестных можно задать произвольно( исключая( может быть( нуль( в зависимости от того( какой параметр взять в качестве произвольного( Положим( например( a2=((0( тогда a1=1-(( b1=b2=1/2( и соотношения 1(9( 1(10( 1(11 сведутся к ym 1=ym h O(h3) 1(12 Это наиболее общая форма записи метода Рунге-Кутта второго порядка( При (=1/2 мы получаем исправленный метод Эйлера( при (=1 получаем модификационный метод Эйлера( Для всех (( отличных от нуля( ошибка ограничения равна e =kh3 1(13Методы Рунге-Кутта третьего и четвертого порядков можно вывести совершенно аналогично тому( как это делалось при выводе методов первого и второго порядков( Мы не будем воспроизводить выкладки( а ограничимся тем( что приведем формулы( описывающие метод четвертого порядка( один из самых употребляемых методов интегрирования дифференциальных уравнений( Этот классический метод Рунге-Кутта описывается системой следующих пяти соотношений ym 1=ym h/6(R1 2R2 2R3 R4) 1(14 где R1=f(xm(ym)( 1(15 R2=f(xm h/2(ym hR1/2)( 1(16 R3=f(xm h/2(ym hR2/2)( 1(17 R4=f(xm h/2(ym hR3/2).
1(18Ошибка ограничения для этого метода равна e =kh5 так что формулы 1(14-1(18 описывают метод четвертого порядка( Заметим( что при использовании этого метода функцию необходимо вычислять четыре раза( 3. Выбор метода реализации программы Исходя из вышеизложенного, для решения систем дифференциальных уравнений мы выбираем наиболее точный метод решения – метод Рунге-Кутта 4 порядка, один из самых употребляемых методов интегрирования дифференциальных уравнений( - этот метод является одноступенчатым и одношаговым - требует информацию только об одной точке - имеет небольшую погрешность - значение функции рассчитывается при каждом шаге 4. Блок-схема программмы Основная программаПроцедура I I Вход f1,C f1,k1,k2,k3,k4 f1,X ,Xk,dp, ,eps,p выход 5. ПрограммаPROGRAM smi h 04; USES cr ; VAR i, :i eger; sum,k1,k2,k3,k4,p,dp,eps,X ,Xk,X,dX:real; rSR,C,dC,r1,r2,r3,r4,cPR:array of real; f1,f2: ex ; PROCEDURE Difur; BEGI dC:=C k2; {dcC} E D; PROCEDURE RK 4; BEGI Difur; FOR i:=1 O DO BEGI r1 (dX/2); E D; Difur; FOR i:=1 O DO BEGI r2 (dX/2); E D; Difur; FOR i:=1 O DO BEGI r3 dX; E D; Difur; FOR i:=1 O DO r4) (r2))/6; E D; PROCEDURE S ROKA; BEGI WRI E(f2,' ',x:4:1,' ',c:7:3,' '); WRI E(f2,sum:3:0,' ',dc:7:3,' '); WRI EL (f2); E D; PROCEDURE RU ; BEGI WRI E('S ep 3: Calcula i g da a a d wri i g resul s o file : ou .rez'); X:=X ; dX:=0.05; REPEA IF (ABS(x-p)Xk); WRI EL (' - do e.'); E D; PROCEDURE I I ; BEGI ClrScr; WRI EL ('Smi h-04: v1.0 (c) 1998 by Mike Smi h smi h01@home.bar.ru '); WRI EL ; WRI EL ; WRI E('S ep 1: Read da a from file : i .da '); ASSIG (f1,'i .da '); RESE (f1); READL (f1,C); READL (f1,k1,k2,k3,k4); READL (f1,X ,Xk,dp, ,eps,p); WRI EL (' - do e.'); ASSIG (f2,'ou .rez'); REWRI E(f2); WRI E('S ep 2: Wri e header o file : ou .rez'); WRI EL (f2,' ,c Ca,% Cb,% Cc,% SUM dCa dCb dCc '); WRI EL (' - do e.'); E D; PROCEDURE DO E; BEGI WRI EL ('S ep 4: Close all files a d exi i g.'); CLOSE(f2); WRI EL ; E D;BEGI I I ; RU ; DO E;E D. 6. Идентификация переменных Таблица 1 7. Результаты расчета Таблица 2 8. Обсуждение результатов расчета. В результате расчета кинетической схемы процесса на языке Паскаль методом Рунге-Кутты, были получены результаты зависимости изменения концентрации реагирующих веществ во времени. Исходя из полученных результатов, можно сделать вывод, что расчет произведен верно, так как, исходя из полученных значений скоростей реакций можно сделать вывод, что соблюдается баланс скоростей химической реакции. Рассмотрим процесс подробнее. Вещество А на протяжении всего процесса расходуется на образование веществ В и С. Концентрации вещества А в начальный момент времени расходуется быстрее, чем концентрации его же в конце процесса. Это обусловлено тем, что скорость химической реакции зависит от концентрации реагирующего вещества. Производная имеет знак «минус». Это говорит о том, что вещество расходуется. Следовательно, чем выше концентрация вещества, вступающего в процесс, тем выше скорость его реагирования с другими веществами. Вещества В и С образуются пропорционально, так как, исходя из кинетической схемы процесса и значений констант скоростей химической реакции, видно, что образование этих веществ и расходование этих веществ, одинаково.
С помощью параметров и опций можно задать множество возможностей для наглядной визуализации довольно сложных решений систем дифференциальных уравнений с частными производными. Следует отметить, что неправильное задание параметров ведет просто к выводу функции в строке вывода без построения графиков и нередко без сообщений об ошибках. Поэтому полезно внимательно просмотреть примеры применения этой функции — как приведенные ниже, так и в справке. 7.8.5. Примеры применения функции PDEplot Рисунок 7.28 демонстрирует применение функции PDEplot. Этот пример из справки показывает, насколько необычным может быть решение даже простой системы дифференциальных уравнений в частных производных. Рис. 7 28. Пример применения функции PDEplot В данном случае решение представлено трехмерной фигурой весьма нерегулярного вида. Другой пример использования функции PDEplot показан на рис. 7.29. Он иллюстрирует комбинированное построение графиков решения разного типа с применением функциональной закраски, реализуемой по заданной формуле с помощью опции initcolor. Рис. 7.29
3. Решение дифференциальных уравнений 1 порядка методом Эйлера
4. Поиск решений системы линейных уравнений методом Гаусса
9. Решение систем дифференциальных уравнений при помощи неявной схемы Адамса 3-го порядка
10. Методы оптимизации при решении уравнений
11. Методы и приемы решения задач
12. Сущность и методы принятия управленческих решений
14. Волновое уравнение не имеет единственного решения
15. Методы поиска технических решений
17. Методы планирования управленческих решений
18. Методы принятия управленческих решений для конкретной проблемы
19. Экспертные методы оценки управленческого решения
20. Методы путевого анализа и их применение к системам одновременных уравнений
21. Системы линейных уравнений
25. Алгебра матриц. Системы линейных уравнений
26. Системы линейных уравнений
27. Методы воспроизведения заболеваний органов системы пищеварения
30. Правовая система России во 2-й половине XlX - начале ХХ вв. Судебная реформа
31. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ПЯТИТОЧЕЧНЫМ МЕТОДОМ АДАМСА – БАШФОРТА
33. Методы и алгоритмы компьютерного решения дифференциальных уравнений
34. Аналитические свойства решений системы двух дифференциальных уравнений третьего порядка
36. Решение систем линейных алгебраических уравнений методом Гаусса и Зейделя
37. Построение решения задачи Гурса для телеграфного уравнения методом Римана
41. Методы решения уравнений в странах древнего мира
42. Приближённые методы решения алгебраического уравнения
43. Решение задачи Дирихле для уравнения Лапласа методом сеток
44. Существование решения дифференциального уравнения и последовательные приближения
45. Метод решения уравнений Ньютона - Рафсона
46. Решение дифференциального уравнения с последующей аппроксимацией
47. Методы решения уравнений, содержащих параметр
48. Метод касательных решения нелинейных уравнений
49. Итерационные методы решения нелинейных уравнений
51. Решение обыкновенных дифференциальных уравнений
52. Решение систем нелинейных алгебраических уравнений методом Ньютона
53. Численные методы решения систем линейных уравнений
57. Решение дифференциальных уравнений
58. Решение систем дифференциальных уравнений
59. Методы решения алгебраических уравнений
61. Прямые методы решения систем линейных алгебраических уравнений
62. Решение транспортной задачи методом потенциалов
64. Билеты, решения и методичка по Информатике (2.0)
65. Краткие сведения о электронных таблицах. Решение уравнения
67. Решение задач - методы спуска
68. Метод последовательных уступок (Теория принятия решений)
69. Методы решения систем линейных неравенств
73. Проект создания системы поддержки принятия решений оперативно-дежурной службы милиции
74. Определение точного коэффициента электропроводности из точного решения кинетического уравнения
75. Методология и методы принятия решения
76. Модели и методы принятия решений
78. Методология и методы принятия решения
80. Системы поддержки и принятия решений
81. Применение графиков в решении уравнений
82. Методы решения некорректно поставленных задач
83. Решение задач линейной оптимизации симплекс – методом
84. Способы решения систем линейных уравнений
85. Алгебраическое и графическое решение уравнений, содержащих модули
90. Эвристические методы решения творческих задач
92. Кинезиология как Метод решения психологических проблем
93. Решение задач транспортного типа методом потенциалов
94. Педагогическая система обучения принятию управленческого решения в вузе
96. Применение новейших экономико-математических методов для решения задач
98. Решение геоэкологических проблем с помощью нестандартных геофизических методов
99. Исследование системы относительно глобальной геополитики и решения конфликта english