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

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

Рубрика Программирование, компьютеры и кибернетика
Вид курсовая работа
Язык русский
Дата добавления 06.01.2013
Размер файла 929,7 K

Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже

Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

РЕФЕРАТ

Курсовая работа: 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, протекающих по следующей схеме:

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

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

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

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=1/2;

b1=b2=1

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

a1=0 a2=1

b1=b2=1/2

Формулы 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=1/2

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

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

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

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

Это наиболее общая форма записи метода Рунге-Кутта второго порядка При =1/2 мы получаем исправленный метод Эйлера при =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_kutt_4;

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_kutt_4;

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

Размещено на Allbest.ru


Подобные документы

Работы в архивах красиво оформлены согласно требованиям ВУЗов и содержат рисунки, диаграммы, формулы и т.д.
PPT, PPTX и PDF-файлы представлены только в архивах.
Рекомендуем скачать работу.