Разработка программы на языке Turbo Pascal 7.0 для решения дифференциальных уравнений
РЕФЕРАТ
Курсовая работа: 25 страниц, 4 рисунка, 2 таблицы, 2 графика, 3 источника
Цель работы: разработать программу на языке Turbo Pascal 7.0 для решения
дифференциальных уравнений.
В курсовой работе приведено описание методов Рунге-Кутта первого, второго,
третьего и четвертого порядков, приведена точность нахождения концентраций
компонентов по рассматриваемым методам, блок-схема алгоритма, разработана
программа расчета на на языке Turbo Pascal 7.0 с результатами и с тестовым
вариантом расчета, проведен анализ полученных результатов.
МЕТОД РУНГЕ-КУТТА, РЕАКЦИЯ, БЛОК-СХЕМА, АЛГОРИТМ, ПРОИЗВОДНАЯ, ТОЧНОСТЬ,
ПРОГРАММА, РАСЧЕТ, ГРАФИК.
СОДЕРЖАНИЕ
ВВЕДЕНИЕ
. ПОСТАНОВКА
ЗАДАЧИ
. СУЩНОСТЬ
МЕТОДОВ РУНГЕ-КУТТА
.
ИДЕНТИФИКАЦИЯ ПЕРЕМЕННЫХ
4. ПРОГРАММА
РАСЧЕТА НА ЯЗЫКЕ TURBO PASCAL
РЕЗУЛЬТАТЫ
РАСЧЕТА
. АНАЛИЗ
РЕЗУЛЬТАТОВ
ЗАКЛЮЧЕНИЕ
ПЕРЕЧЕНЬ
ССЫЛОК
ВВЕДЕНИЕ
Многие процессы химической технологии описываются системой
дифференциальных уравнений - начиная от кинетических исследований и заканчивая
химическими технологическими процессами. В основу математических способов описания процессов положены
система дифференциальных уравнений и система линейных алгебраических уравнений. Эти уравнения описывают материальные
и тепловые балансы объектов химической технологии, а так же структуры потоков технических веществ в этих
аппаратах.
Для получения,
распределения технологических параметров во времени и в пространстве (в
пределах объекта),
необходимо произвести СДУ методом, который дал бы высокую точность решения при минимальных затратах
времени на решение,
потому что ЭВМ должна работать в режиме реального времени и успевать за ходом
технологического процесса. Если время на решение задачи большое, то управляющее воздействие, выработанное на ЭВМ может привести к
отрицательным воздействиям. Методов решения существует очень много. В данной работе будет изучена кинетическая схема
протекания химических реакций, которая будет переведена на математический язык
- систему дифференциальных уравнений первого порядка, которая будет решаться в численном
виде методом Рунге-Кутта четвертого порядка.
программа дифференциальные уравнение кинетический
1. ПОСТАНОВКА ЗАДАЧИ
Необходимо найти распределение концентраций компонентов во времени, а
также изменение скорости химических реакций до веществ А,B,C,E,D, протекающих по следующей схеме:
Так
как коэффициенты 
являются константами, то можно уравнение записать в
следующем виде.
Для преобразования данных дифференциальных уравнений
для использования их в расчетах тепловых и кинетических схем методами
Рунге-Кутта необходимо подставлять вместо производных значений концентраций,
значения концентраций данных в начале процесса. Это обусловлено тем, что метод
Рунге-Кутта четвертого порядка, который будет использован для расчета
кинетической схемы процесса. Так как этот метод требует сведений только об
одной точке и значений функции.
2. СУЩНОСТЬ МЕТОДОВ РУНГЕ-КУТТА
Разбор и рассмотрение методов, применяемых на практике для решения дифференциальных
уравнений, мы начнем с их
широкой категории,
известной под общим названием методов Рунге-Кутта.
Методы Рунге-Кутта обладают следующими свойствами:
a) Эти методы являются одноступенчатыми: чтобы найти уm+1, нужна информация о предыдущей точке
xm,ym.) Они согласуются с рядом Тейлора вплоть до членов порядка hp, где степень р различна для различных методов и
называется порядковым номером или порядком метода.) Они не требуют вычисления производных от 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
тогда уравнение примет вид
m+1=ym+h*f(xm,ym) 1.1
Ошибка при x=xm+1 показана в виде отрезка е. Очевидно, найденное таким образом приближенное
значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h, так что ошибка ограничения равна et=Кh2
Заметим, что хотя точка
на рисунке 2.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. В этой точке снова вычисляется
тангенс, дает прямую L. Наконец, через точку xm,ym мы проводим прямую L, параллельную L. Точка, в которой прямая L пересечется с ординатой, восстановленной из x=xm+1=xm+h, и будет искомой точкой xm+1,ym+1.
Тангенс угла наклона прямой L и прямой L равен
Ф(xm,ym,h)=?[f(xm,ym)+f(xm+h,ym+y¢mh)] 1.2
где y¢m=f(xm,ym) 1.3
Уравнение линии L при этом записывается в виде
=ym+(x-xm)Ф(xm,ym,h),
так что
ym+1=ym+hФ(xm,ym,h). 1.4
Соотношения 1.2, 1.3, 1.4 описывают исправленный метод Эйлера.
Чтобы выяснить,
насколько хорошо этот метод согласуется с разложением в ряд Тейлора, вспомним, что разложение в ряд функции f(x,y) можно записать следующим образом:
(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, получаем
(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 и сравним с рядом Тейлора
m+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
Прямая с таким наклоном, проходящая через Р,
обозначена через L*. Вслед за тем, мы
проводим через точку xm,ym
прямую параллельную L*, и обозначаем ее через L0. Пересечение этой прямой с ординатой x=xm+h и даст
искомую точку xm+1,ym+1.
Уравнение прямой можно записать в виде
y=ym+(x-xm)Ф(xm,ym,h),
где Ф задается формулой 1.6. Поэтому
m+1=ym+hФ(xm,ym,h) 1.8
Соотношения 1.6, 1.7, 1.8 описывают так называемый
модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго
порядка. Обобщим оба метода. Заметим, что оба метода описываются формулами
вида
m+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
В частности,
для исправленного метода Эйлера1=a2=1/2;1=b2=1.
Формулы
1.9, 1.10, 1.11 описывают некоторый метод типа Рунге-Кутта. Посмотрим, какого порядка метод можно рассчитывать получить в
лучшем случае и каковы допустимые значения параметров a1, a2, b1 и b2 .
Чтобы
получить соответствие ряду Тейлора вплоть до членов степени h, в общем случае достаточно одного параметра. Чтобы получить согласование вплоть до членов степени h2, потребуется
еще два параметра, так как необходимо учитывать члены h2fx и h2ffy. Так как у нас имеется всего четыре параметра, три из которых потребуются для создания согласования с рядом Тейлора
вплоть до членов порядка h2, то самое лучшее, на что здесь можно рассчитывать - это метод второго
порядка.
В
разложении f(x,y) в ряд 1.5 в окрестности
точки xm,ym положим
=xm+b1h,=ym+b2hf.
Тогда
(xm+b1h,ym+b2hf)=f+b1hfx+b2hffy+O(h2),
где
функция и производные в правой части равенства вычислены в точке xm,ym.
Тогда
1.9 можно переписать в виде
+1=ym+h[a1f+a2f+h(a2b1fx+a2b2ffy)]+O(h3).
Сравнив
эту формулу с разложением в ряд Тейлора, можно переписать в
виде
m+1=ym+h[a1f+a2f+h(a2b1fx+a2b2ffy)]+O(h3).
Если
потребовать совпадения членов hf, то a1+a2=1.
Сравнивая
члены, содержащие h2fx, получаем a2b1=1/2.
Сравнивая
члены, содержащие h2ffy, получаем a2b2=1/2.
Так
как мы пришли к трем уравнениям для определения четырех неизвестных, то одно из этих неизвестных можно задать произвольно, исключая, может быть, нуль, в зависимости от того, какой параметр взять в качестве
произвольного.
Положим, например, a2=w¹0. тогда a1=1-w, b1=b2=1/2w и соотношения 1.9, 1.10, 1.11 сведутся к
m+1=ym+h[(1-w)f(xm,ym)+wf(xm+h/2w,ym+h/2wf(xm,ym))]+O(h3) 1.12
Это
наиболее общая форма записи метода Рунге-Кутта второго порядка. При w=1/2 мы получаем исправленный метод Эйлера, при w=1 получаем модификационный метод Эйлера. Для всех w, отличных от нуля, ошибка ограничения
равна
t=kh3 1.13
Методы
Рунге-Кутта третьего и четвертого порядков можно вывести совершенно аналогично
тому, как это делалось при выводе методов первого и второго
порядков. Мы не будем воспроизводить выкладки, а ограничимся тем, что приведем формулы, описывающие метод четвертого порядка, один из
самых употребляемых методов интегрирования дифференциальных уравнений. Этот классический метод Рунге-Кутта описывается системой следующих
пяти соотношений
m+1=ym+h/6(R1+2R2+2R3+R4)
1.14
где
R1=f(xm,ym), 1.152=f(xm+h/2,ym+hR1/2), 1.163=f(xm+h/2,ym+hR2/2), 1.174=f(xm+h/2,ym+hR3/2).
1.18
Ошибка ограничения для этого метода равна et=kh5 так, что формулы 1.14-1.18 описывают метод четвертого порядка. Заметим, что при использовании этого метода
функцию необходимо вычислять четыре раза.
Исходя из вышеизложенного, для решения систем дифференциальных уравнений
мы выбираем наиболее точный метод решения - метод Рунге-Кутта четвертого
порядка, один из самых употребляемых методов интегрирования дифференциальных
уравнений.
Этот метод является одноступенчатым и одношаговым, требует информацию
только об одной точке имеет небольшую погрешность значение функции рассчитывается
при каждом шаге.
Рисунок 2.4 Метод Рунге-Кутта четвертого порядка
3. ИДЕНТИФИКАЦИЯ ПЕРЕМЕННЫХ
Таблица 3.1
. ПРОГРАММА РАСЧЕТА НА ЯЗЫКЕ TURBO PASCAL 7.0
program Kurs;
n=5;c,dc,cpr:array[1..n] of real;:array [1..5] of
real;,p,dp,tn,tk,t,eps:real;:integer;:text;DF;[1]:=-k[1]*c[1];[2]:=k[1]*c[1]-(k[2]+k[5])*c[2]+k[4]*c[3];[3]:=k[2]*c[2]-(k[4]+k[3])*c[3];[4]:=k[3]*c[3];[5]:=k[5]*c[2];;
Runge_kutt_4;r:array[1..5,1..n] of real;;i:= 1 to n do
begin[1,i]:=dc[i];[i]:=cpr[i]+r[1,i]*h/2;;:=t+h/2;;i:= 1 to n do
begin[2,i]:=dc[i];[i]:=cpr[i]+r[3,i]*h/2;;;i:= 1 to n do
begin[3,i]:=dc[i];[i]:=cpr[i]+r[3,i]*h;;:=t+h/2;;i:= 1 to n do
begin[4,i]:=dc[i];[5,i]:=1/6*(r[1,i]+2*r[2,i]+2*r[3,i]+r[4,i]);[i]:=cpr[i]+r[5,i]*h/2;;;
outdata;(f,t:1:4,#9,c[1]:1:4,#9,dc[1]:1:4,#9,c[2]:1:4,#9,dc[2]:1:4,#9,c[3]:1:4,#9,dc[3]:1:4,#9,[4]:1:4,#9,dc[4]:1:4,#9,c[5]:1:4,#9,dc[5]:1:4,#9,c[1]+c[2]+c[3]+c[4]+c[5]:1:4);;
(f,'in.txt');(f);(f,c[1],c[2],c[3],c[4],c[5]);(f,k[1],k[2],k[3],k[4],k[5]);(f,h,dp,tn,tk,eps);(f);
(f,'out.txt');(f);
:=tn;:=tn;('-------------Begin data--------------');('| c1 = ',c[1]:10:4,
' | k1 = ',k[1]:10:4, ' | ');('| c2 = ',c[2]:10:4, ' | k2 = ',k[2]:10:4, ' |
');('| c3 = ',c[3]:10:4, ' | k3 = ',k[3]:10:4, ' | ');('| c4 = ',c[4]:10:4, ' |
k4 = ',k[4]:10:4, ' | ');('| c5 = ',c[5]:10:4, ' | k5 = ',k[5]:10:4, ' |
');('> sum ',c[1]+c[2]+c[3]+c[4]+c[5]:10:4,' |');;('| h = ', h:10:4, ' | dp
= ', dp:10:4, ' | ');('| tn = ', tn:10:4, ' | tk = ', tk:10:4, ' | ');('|eps =
', eps:10:4, ' | |');('------------End data-----------------');;
(f,'t',#9,'c[1]',#9,'dc[1]',#9,'c[2]',#9,'dc[2]',#9,'c[3]',#9,'dc[3]',#9,
'c[4]',#9,'dc[4]',#9,'c[5]',#9,'dc[5]',#9,'sum');
abs(t-p)<eps then begin;p=20 then;;:=p+dp;;i:= 1 to n do
cpr[i]:=c[i];_kutt_4;t>tk+eps;(f);;.
Пример файла исходных данных in.txt
1 1 0 0
.3 0.2 0.05 0.1 0.05
.01.5 0 20 0.001
. РЕЗУЛЬТАТЫ РАСЧЕТА
Таблица 5.1
6. АНАЛИЗ РЕЗУЛЬТАТОВ РАСЧЕТА
Вещество А на протяжении всего процесса расходуется на образование
веществ В, С, E, D. Концентрации вещества А в начальный момент времени
расходуется быстрее, чем концентрации его же в конце процесса. Это обусловлено
тем, что скорость химической реакции зависит от концентрации реагирующего
вещества. Производная имеет знак «минус». Это говорит о том, что вещество
расходуется. Следовательно, чем выше концентрация вещества, вступающего в процесс,
тем выше скорость его реагирования с другими веществами. Вещество В образуется
из вещества A и С и расходуется на производство веществ D и С. Поскольку
концентрация вещества А большая в начале процесса и расходуется со временем, то
и концентрация вещества В сначала быстро возрастает, а затем скорость
уменьшается и меняет знак с «плюс» на «минус», в результате чего концентрация
вещества начинает медленно уменьшаться.
Аналогичная ситуация происходит с веществом С, но пик его концентрации
приходится наболее дальний срок (уже после 20 сек.). Концентрации веществ D и Е
на протяжении всего времени реакции имеют положительные производные, в
следствии чего можно и концентрации данных веществ ростут со временем. Причем,
поскольку концентрация вещества В, в начале процесса, на много выше, чем
концентрация вещества С, то и концентрация вещества D в процессе реакции будет
всегда выше, чем концентрация вещества E. Реакция протекает в течении
длительного времени и не заканчивается после 20 сек
График 6.1
График 6.2
ЗАКЛЮЧЕНИЕ
В ходе выполнения работы был произведен расчет системы дифференциальных
уравнений методом Рунге-Кутта четвертого порядка, произведен расчет
кинетической схемы процесса при изотермических условиях при данных значениях
концентраций и констант скоростей. Расчет произведен с малой величиной
погрешности.
В результате выполнения расчета получена зависимость изменения
концентрации вещества во времени. Из расчета следует, что на протяжении всего
процесса вещество А расходовалось на образование остальных. Процесс не достиг
конечного состояния (не достиг равновесия).
ПЕРЕЧЕНЬ ССЫЛОК
1 Мудров А.Е.
Численные методы для ПЭВМ на языках Паскаль, Фортран и Бейсик. МП “Раско”,
Томск, 1991 г.
Самарский
А.А. «Введение в численные методы» М.: «Наука» 1978, 269 стр.
Гулин А.В.,
Самарский А.А. «Численные методы» М.: «Наука» 1989, 432стр.