Разработка программы на языке Turbo Pascal 7.0 для решения дифференциальных уравнений

РЕФЕРАТ Курсовая работа: 25 страниц, 4 рисунка, 2 таблицы, 2 графика, 3 источника Цель работы: разработать программу на языке Turbo Pascal 7.0 для решения дифференциальных уравнений.

В курсовой работе приведено описание методов Рунге-Кутта первого, второго, третьего и четвертого порядков, приведена точность нахождения концентраций компонентов по рассматриваемым методам, блок-схема алгоритма, разработана программа расчета на на языке Turbo Pascal 7.0 с результатами и с тестовым вариантом расчета, проведен анализ полученных результатов.

МЕТОД РУНГЕ-КУТТА, РЕАКЦИЯ, БЛОК-СХЕМА, АЛГОРИТМ, ПРОИЗВОДНАЯ, ТОЧНОСТЬ, ПРОГРАММА, РАСЧЕТ, ГРАФИК.

СОДЕРЖАНИЕ ВВЕДЕНИЕ

1. ПОСТАНОВКА ЗАДАЧИ

2. СУЩНОСТЬ МЕТОДОВ РУНГЕ-КУТТА

3. ИДЕНТИФИКАЦИЯ ПЕРЕМЕННЫХ

4. ПРОГРАММА РАСЧЕТА НА ЯЗЫКЕ TURBO PASCAL

5 РЕЗУЛЬТАТЫ РАСЧЕТА

6. АНАЛИЗ РЕЗУЛЬТАТОВ ЗАКЛЮЧЕНИЕ ПЕРЕЧЕНЬ ССЫЛОК

ВВЕДЕНИЕ

Многие процессы химической технологии описываются системой дифференциальных уравнений — начиная от кинетических исследований и заканчивая химическими технологическими процессами В основу математических способов описания процессов положены система дифференциальных уравнений и система линейных алгебраических уравнений Эти уравнения описывают материальные и тепловые балансы объектов химической технологии, а так же структуры потоков технических веществ в этих аппаратах Для получения распределения технологических параметров во времени и в пространстве (в пределах объекта) необходимо произвести СДУ методом который дал бы высокую точность решения при минимальных затратах времени на решение потому что ЭВМ должна работать в режиме реального времени и успевать за ходом технологического процесса. Если время на решение задачи большое то управляющее воздействие выработанное на ЭВМ может привести к отрицательным воздействиям Методов решения существует очень много В данной работе будет изучена кинетическая схема протекания химических реакций, которая будет переведена на математический язык — систему дифференциальных уравнений первого порядка, которая будет решаться в численном виде методом Рунге-Кутта четвертого порядка.

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

1. ПОСТАНОВКА ЗАДАЧИ

Необходимо найти распределение концентраций компонентов во времени, а также изменение скорости химических реакций до веществ А,B, C, E, D, протекающих по следующей схеме:

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

Для преобразования данных дифференциальных уравнений для использования их в расчетах тепловых и кинетических схем методами Рунге-Кутта необходимо подставлять вместо производных значений концентраций, значения концентраций данных в начале процесса. Это обусловлено тем, что метод Рунге-Кутта четвертого порядка, который будет использован для расчета кинетической схемы процесса. Так как этот метод требует сведений только об одной точке и значений функции.

2. СУЩНОСТЬ МЕТОДОВ РУНГЕ-КУТТА

Разбор и рассмотрение методов применяемых на практике для решения дифференциальных уравнений мы начнем с их широкой категории известной под общим названием методов Рунге-Кутта Методы Рунге-Кутта обладают следующими свойствами:

a) Эти методы являются одноступенчатыми: чтобы найти у_m+1 нужна информация о предыдущей точке xmym

b) Они согласуются с рядом Тейлора вплоть до членов порядка hp где степень р различна для различных методов и называется порядковым номером или порядком метода

c) Они не требуют вычисления производных от f (xy) а требуют вычисления самой функции Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий После этого мы подтвердим полученные результаты аналитически Предположим нам известна точка xmym на искомой кривой Тогда мы можем провести прямую линию с тангенсом угла наклона у_m=f (x_my_m) которая пройдет через точку x_my_m Это построение показано на рис1 где кривая представляет собой точное, но конечно неизвестное решение уравнения, а прямая линия L1 построена так как это только что описано Тогда следующей точкой решения можно считать ту где прямая L₁ пересечет ординату проведенную через точку x=x_m+1=x_m+h

Уравнение прямой L₁ выглядит так: y=y_m+y_m(x-x_m) так как y=f (x_my_m) и кроме того x_m+1=x_m+h тогда уравнение примет вид

y_m+1=y_m+h*f (x_my_m) 11

Ошибка при x=x_m+1 показана в виде отрезка е Очевидно найденное таким образом приближенное значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h так что ошибка ограничения равна e_t=Кh2

Заметим что хотя точка на рисунке 2.1 была показана на кривой в действительности y_m является приближенным значением и не лежит точно на кривой Формула 11 описывает метод Эйлера один из самых старых и широко известных методов численного интегрирования дифференциальных уравнений Отметим что метод Эйлера является одним из методов Рунге-Кутта первого порядка Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: x_my_m и x_m+hy_m+hy_m Последняя точка есть та самая которая в методе Эйлера обозначалась x_m+1y_m+1 Геометрический процесс нахождения точки x_m+1y_m+1 можно проследить по рис2 С помощью метода Эйлера находится точка x_m+hy_m+hy_m лежащая на прямой L₁ В этой точке снова вычисляется тангенс дает прямую Наконец через точку x_my_m мы проводим прямую L параллельную Точка в которой прямая L пересечется с ординатой восстановленной из x=x_m+1=x_m+h и будет искомой точкой x_m+1y_m+1

Тангенс угла наклона прямой и прямой L равен

Ф (x_my_mh)=[f (x_my_m)+f (x_m+hy_m+y_mh)] 12

где y_m=f (x_my_m) 13

Уравнение линии L при этом записывается в виде

y=y_m+(x-x_m)Ф (x_my_mh)

так что

y_m+1=y_m+hФ (x_my_mh) 14

Соотношения 12 13 14 описывают исправленный метод Эйлера Чтобы выяснить насколько хорошо этот метод согласуется с разложением в ряд Тейлора вспомним что разложение в ряд функции f (xy) можно записать следующим образом:

f (xy)=f (x_my_m)+(x-x_m)f/x+(y-y_m)f/x+ 15

где частные производные вычисляются при x=x_m и y=y_m

Подставляя в формулу 15 x=x_m+h и y=y_m+hy_m и используя выражение 13 для y_m получаем

f (x_m+hy_m+hy_m)=f+hf_x+hff_y+O (h2)

где снова функция f и ее производные вычисляются в точке x_my_m Подставляя результат в 12 и производя необходимые преобразования получаем Ф (x_my_mh)=f+h/2(f_x+ff_y)+O (h2)

Подставим полученное выражение в 14 и сравним с рядом Тейлора

y_m+1=y_m+hf+h2/2(f_x+ff_y)+O (h3)

Как видим исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h2 являясь таким образом методом Рунге-Кутты второго порядка Рассмотрим модификационный метод Эйлера Рассмотрим рис3 где первоначальное построение сделано так же как и на рис2 Но на этот раз мы берем точку лежащую на пересечении этой прямой и ординатой x=x+h/2 На рисунке эта точка образована через Р, а ее ордината равна y=y_m+(h/2)y_m Вычислим тангенс угла наклона касательной в этой точке Ф (x_my_mh)=f+(x_m+h/2y_m+h/2*y_m) 16

где y_m=f (x_my_m) 17

Прямая с таким наклоном проходящая через Р обозначена через * Вслед за тем мы проводим через точку xmym прямую параллельную * и обозначаем ее через L₀ Пересечение этой прямой с ординатой x=x_m+h и даст искомую точку x_m+1y_m+1 Уравнение прямой можно записать в виде

y=y_m+(x-x_m)Ф (x_my_mh)

где Ф задается формулой 16 Поэтому

y_m+1=y_m+hФ (x_my_mh) 18

Соотношения 16 17 18 описывают так называемый модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго порядка Обобщим оба метода Заметим что оба метода описываются формулами вида

y_m+1=y_m+hФ (x_my_mh) 19

и в обоих случаях Ф имеет вид Ф (x_my_mh)=a₁f (x_my_m)+a₂f (x_m+b₁hy_m+b₂hy_m) 110

где y_m=f (x_my_m) 111

В частности для исправленного метода Эйлера

a₁=a₂=½;

b₁=b₂=1

В то время как для модификационного метода Эйлера

a₁=0 a₂=1

b₁=b₂=½

Формулы 19 110 111 описывают некоторый метод типа Рунге-Кутта Посмотрим какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a₁ a₂ b₁ и b₂

Чтобы получить соответствие ряду Тейлора вплоть до членов степени h в общем случае достаточно одного параметра Чтобы получить согласование вплоть до членов степени h2 потребуется еще два параметра так как необходимо учитывать члены h2f_x и h2ff_y Так как у нас имеется всего четыре параметра три из которых потребуются для создания согласования с рядом Тейлора вплоть до членов порядка h2 то самое лучшее на что здесь можно рассчитывать — это метод второго порядка В разложении f (xy) в ряд 15 в окрестности точки x_my_m положим

x=x_m+b₁h

y=y_m+b₂hf

Тогда

f (x_m+b₁hy_m+b₂hf)=f+b₁hf_x+b₂hff_y+O (h2)

где функция и производные в правой части равенства вычислены в точке x_my_m

Тогда 19 можно переписать в виде

ym+1=ym+h[a1f+a2f+h (a2b1fx+a2b2ffy)]+O (h3)

Сравнив эту формулу с разложением в ряд Тейлора можно переписать в виде

y_m+1=y_m+h[a₁f+a₂f+h (a₂b₁f_x+a₂b₂ff_y)]+O (h3)

Если потребовать совпадения членов hf то a₁+a₂=1

Сравнивая члены содержащие h2f_x получаем a₂b₁=½

Сравнивая члены содержащие h2ff_y получаем a₂b₂=½

Так как мы пришли к трем уравнениям для определения четырех неизвестных то одно из этих неизвестных можно задать произвольно исключая может быть нуль в зависимости от того какой параметр взять в качестве произвольного

Положим например a₂=0 тогда a₁=1- b₁=b₂=½ и соотношения 19 110 111 сведутся к

y_m+1=y_m+h[(1-)f (x_my_m)+f (x_m+h/2y_m+h/2f (x_my_m))]+O (h3) 112

Это наиболее общая форма записи метода Рунге-Кутта второго порядка При =½ мы получаем исправленный метод Эйлера при =1 получаем модификационный метод Эйлера Для всех отличных от нуля ошибка ограничения равна

e_t=kh3 113

Методы Рунге-Кутта третьего и четвертого порядков можно вывести совершенно аналогично тому как это делалось при выводе методов первого и второго порядков Мы не будем воспроизводить выкладки, а ограничимся тем что приведем формулы описывающие метод четвертого порядка один из самых употребляемых методов интегрирования дифференциальных уравнений Этот классический метод Рунге-Кутта описывается системой следующих пяти соотношений

y_m+1=y_m+h/6(R₁+2R₂+2R₃+R₄) 114

где R₁=f (x_my_m) 115

R₂=f (x_m+h/2y_m+hR₁/2) 116

R₃=f (x_m+h/2y_m+hR₂/2) 117

R₄=f (x_m+h/2y_m+hR₃/2). 118

Ошибка ограничения для этого метода равна e_t=kh5 так, что формулы 114−118 описывают метод четвертого порядка Заметим что при использовании этого метода функцию необходимо вычислять четыре раза Исходя из вышеизложенного, для решения систем дифференциальных уравнений мы выбираем наиболее точный метод решения — метод Рунге-Кутта четвертого порядка, один из самых употребляемых методов интегрирования дифференциальных уравнений

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

Рисунок 2.4 Метод Рунге-Кутта четвертого порядка

3. ИДЕНТИФИКАЦИЯ ПЕРЕМЕННЫХ

Таблица 3.1

4. ПРОГРАММА РАСЧЕТА НА ЯЗЫКЕ TURBO PASCAL 7.0

program Kurs;

const n=5;

var c, dc, cpr: array[1.n] of real;

k:array [1.5] of real;

h, p, dp, tn, tk, t, eps:real;

i:integer;

f:text;

procedure DF;

begin

dc[1]: =-k[1]*c[1];

dc[2]:=k[1]*c[1]-(k[2]+k[5])*c[2]+k[4]*c[3];

dc[3]:=k[2]*c[2]-(k[4]+k[3])*c[3];

dc[4]:=k[3]*c[3];

dc[5]:=k[5]*c[2];

end;

Procedure Runge_kutt₄;

var r: array[1.5,1.n] of real;

begin

DF;

for i:= 1 to n do begin

r[1,i]: =dc[i];

c[i]:=cpr[i]+r[1,i]*h/2;

end;

t:=t+h/2;

DF;

for i:= 1 to n do begin

r[2,i]: =dc[i];

c[i]:=cpr[i]+r[3,i]*h/2;

end;

DF;

for i:= 1 to n do begin

r[3,i]: =dc[i];

c[i]:=cpr[i]+r[3,i]*h;

end;

t:=t+h/2;

DF;

for i:= 1 to n do begin

r[4,i]: =dc[i];

r[5,i]:=1/6*(r[1,i]+2*r[2,i]+2*r[3,i]+r[4,i]);

c[i]:=cpr[i]+r[5,i]*h/2;

end;

end;

procedure outdata;

begin

writeln (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,

c[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);

end;

begin

Assign (f,'in.txt');

Reset (f);

readln (f, c[1], c[2], c[3], c[4], c[5]);

readln (f, k[1], k[2], k[3], k[4], k[5]);

readln (f, h, dp, tn, tk, eps);

Close (f);

Assign (f,'out.txt');

Rewrite (f);

t:=tn;

p:=tn;

writeln ('——————-Begin data———————');

writeln ('| c1 = ', c[1]: 10:4, ' | k1 = ', k[1]: 10:4, ' | ');

writeln ('| c2 = ', c[2]: 10:4, ' | k2 = ', k[2]: 10:4, ' | ');

writeln ('| c3 = ', c[3]: 10:4, ' | k3 = ', k[3]: 10:4, ' | ');

writeln ('| c4 = ', c[4]: 10:4, ' | k4 = ', k[4]: 10:4, ' | ');

writeln ('| c5 = ', c[5]: 10:4, ' | k5 = ', k[5]: 10:4, ' | ');

writeln ('> sum ', c[1]+c[2]+c[3]+c[4]+c[5]: 10:4,' |');

writeln;

writeln ('| h = ', h:10:4, ' | dp = ', dp:10:4, ' | ');

writeln ('| tn = ', tn:10:4, ' | tk = ', tk:10:4, ' | ');

writeln ('|eps = ', eps:10:4, ' | |');

writeln ('——————End data————————-');

writeln;

writeln (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');

repeat

if abs (t-p)

DF;

if p=20 then

outdata;

outdata;

P:=p+dp;

end;

for i:= 1 to n do cpr[i]: =c[i];

Runge_kutt₄;

until t>tk+eps;

close (f);

ReadKey;

end.

Пример файла исходных данных in. txt

98 1 1 0 0

0.3 0.2 0.05 0.1 0.05

0.01.5 0 20 0.001

5. РЕЗУЛЬТАТЫ РАСЧЕТА Таблица 5.1

6. АНАЛИЗ РЕЗУЛЬТАТОВ РАСЧЕТА Вещество, А на протяжении всего процесса расходуется на образование веществ В, С, E, D. Концентрации вещества, А в начальный момент времени расходуется быстрее, чем концентрации его же в конце процесса. Это обусловлено тем, что скорость химической реакции зависит от концентрации реагирующего вещества. Производная имеет знак «минус». Это говорит о том, что вещество расходуется. Следовательно, чем выше концентрация вещества, вступающего в процесс, тем выше скорость его реагирования с другими веществами. Вещество В образуется из вещества A и С и расходуется на производство веществ D и С. Поскольку концентрация вещества, А большая в начале процесса и расходуется со временем, то и концентрация вещества В сначала быстро возрастает, а затем скорость уменьшается и меняет знак с «плюс» на «минус», в результате чего концентрация вещества начинает медленно уменьшаться.

Аналогичная ситуация происходит с веществом С, но пик его концентрации приходится наболее дальний срок (уже после 20 сек.). Концентрации веществ D и Е на протяжении всего времени реакции имеют положительные производные, в следствии чего можно и концентрации данных веществ ростут со временем. Причем, поскольку концентрация вещества В, в начале процесса, на много выше, чем концентрация вещества С, то и концентрация вещества D в процессе реакции будет всегда выше, чем концентрация вещества E. Реакция протекает в течении длительного времени и не заканчивается после 20 сек

График 6.1

График 6.2

ЗАКЛЮЧЕНИЕ

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

В результате выполнения расчета получена зависимость изменения концентрации вещества во времени. Из расчета следует, что на протяжении всего процесса вещество, А расходовалось на образование остальных. Процесс не достиг конечного состояния (не достиг равновесия).

ПЕРЕЧЕНЬ ССЫЛОК

1 Мудров А. Е. Численные методы для ПЭВМ на языках Паскаль, Фортран и Бейсик. МП «Раско», Томск, 1991 г.

2 Самарский А. А. «Введение в численные методы» М.: «Наука» 1978, 269 стр.

3 Гулин А. В., Самарский А. А. «Численные методы» М.: «Наука» 1989, 432стр.