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

Размещено на 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(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₂=1/2;

b₁=b₂=1

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

a₁=0 a₂=1

b₁=b₂=1/2



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

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

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

Положим например a₂=0 тогда a₁=1- b₁=b₂=1/2 и соотношения 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/2 мы получаем исправленный метод Эйлера при =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_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