Разработка программы на языке 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 (xmym) которая пройдет через точку xmym Это построение показано на рис1 где кривая представляет собой точное, но конечно неизвестное решение уравнения, а прямая линия L1 построена так как это только что описано

Тогда следующей точкой решения можно считать ту где прямая L1 пересечет ординату проведенную через точку x=xm+1=xm+h

Уравнение прямой L1 выглядит так: y=ym+ym(x-xm) так как y=f (xmym) и кроме того xm+1=xm+h тогда уравнение примет вид

ym+1=ym+h*f (xmym) 11

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

Заметим что хотя точка на рисунке 2.1 была показана на кривой в действительности ym является приближенным значением и не лежит точно на кривой

Формула 11 описывает метод Эйлера один из самых старых и широко известных методов численного интегрирования дифференциальных уравнений Отметим что метод Эйлера является одним из методов Рунге-Кутта первого порядка

Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: xmym и xm+hym+hym Последняя точка есть та самая которая в методе Эйлера обозначалась xm+1ym+1 Геометрический процесс нахождения точки xm+1ym+1 можно проследить по рис2 С помощью метода Эйлера находится точка xm+hym+hym лежащая на прямой L1 В этой точке снова вычисляется тангенс дает прямую Наконец через точку xmym мы проводим прямую L параллельную Точка в которой прямая L пересечется с ординатой восстановленной из x=xm+1=xm+h и будет искомой точкой xm+1ym+1

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

Ф (xmymh)=[f (xmym)+f (xm+hym+ymh)] 12

где ym=f (xmym) 13

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

y=ym+(x-xm)Ф (xmymh)

так что

ym+1=ym+hФ (xmymh) 14

Соотношения 12 13 14 описывают исправленный метод Эйлера

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

f (xy)=f (xmym)+(x-xm)f/x+(y-ym)f/x+ 15

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

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

f (xm+hym+hym)=f+hfx+hffy+O (h2)

где снова функция f и ее производные вычисляются в точке xmym Подставляя результат в 12 и производя необходимые преобразования получаем

Ф (xmymh)=f+h/2(fx+ffy)+O (h2)

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

ym+1=ym+hf+h2/2(fx+ffy)+O (h3)

Как видим исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h2 являясь таким образом методом Рунге-Кутты второго порядка

Рассмотрим модификационный метод Эйлера Рассмотрим рис3 где первоначальное построение сделано так же как и на рис2 Но на этот раз мы берем точку лежащую на пересечении этой прямой и ординатой x=x+h/2 На рисунке эта точка образована через Р, а ее ордината равна y=ym+(h/2)ym Вычислим тангенс угла наклона касательной в этой точке

Ф (xmymh)=f+(xm+h/2ym+h/2*ym) 16

где ym=f (xmym) 17

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

y=ym+(x-xm)Ф (xmymh)

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

ym+1=ym+hФ (xmymh) 18

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

ym+1=ym+hФ (xmymh) 19

и в обоих случаях Ф имеет вид

Ф (xmymh)=a1f (xmym)+a2f (xm+b1hym+b2hym) 110

где ym=f (xmym) 111

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

a1=a2=½;

b1=b2=1

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

a1=0 a2=1

b1=b2

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

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

В разложении f (xy) в ряд 15 в окрестности точки xmym положим

x=xm+b1h

y=ym+b2hf

Тогда

f (xm+b1hym+b2hf)=f+b1hfx+b2hffy+O (h2)

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

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

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

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

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

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

Сравнивая члены содержащие h2fx получаем a2b1

Сравнивая члены содержащие h2ffy получаем a2b2

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

Положим например a2=0 тогда a1=1- b1=b2=½ и соотношения 19 110 111 сведутся к

ym+1=ym+h[(1-)f (xmym)+f (xm+h/2ym+h/2f (xmym))]+O (h3) 112

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

et=kh3 113

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

ym+1=ym+h/6(R1+2R2+2R3+R4) 114

где R1=f (xmym) 115

R2=f (xm+h/2ym+hR1/2) 116

R3=f (xm+h/2ym+hR2/2) 117

R4=f (xm+h/2ym+hR3/2). 118

Ошибка ограничения для этого метода равна et=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_kutt4;

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)< eps then begin

DF;

if p=20 then

outdata;

outdata;

P: =p+dp;

end;

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

Runge_kutt4;

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стр.

ПоказатьСвернуть
Заполнить форму текущей работой