Решение систем дифференциальных уравнений методом Рунге-Куты 4 порядка - (курсовая)
p>Соотношения 1. 2, 1. 3, 1. 4 описывают исправленный метод Эйлера. Чтобы выяснить, насколько хорошо этот метод согласуется с разложением в ряд Тейлора, вспомним, что разложение в ряд функции f(x, y) можно записать следующим образом:f(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, получаем f(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 и сравним с рядом Тейлора ym+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. Поэтому
ym+1=ym+hФ(xm, ym, h) 1. 8
Соотношения 1. 6, 1. 7, 1. 8 описывают так называемый модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго порядка. Обобщим оба метода. Заметим, что оба метода описываются формулами вида
ym+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
В частности, для исправленного метода Эйлера
a1=a2=1/2;
b1=b2=1.
В то время как для модификационного метода Эйлера
a1=0, a2=1,
b1=b2=1/2.
Формулы 1. 9, 1. 10, 1. 11 описывают некоторый метод типа Рунге-Кутты. Посмотрим, какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a1, a2, b1 и b2 .
Чтобы получить соответствие ряду Тейлора вплоть до членов степени h, в общем случае достаточно одного параметра. Чтобы получить согласование вплоть до членов степени h2, потребуется еще два параметра, так как необходимо учитывать члены h2fx и h2ffy. Так как у нас имеется всего четыре параметра, три из которых потребуются для создания согласования с рядом Тейлора вплоть до членов порядка h2, то самое лучшее, на что здесь можно рассчитывать - это метод второго порядка.
В разложении f(x, y) в ряд 1. 5 в окрестности точки xm, ym положим x=xm+b1h, y=ym+b2hf.
Тогда f(xm+b1h, ym+b2hf)=f+b1hfx+b2hffy+O(h2), где функция и производные в правой части равенства вычислены в точке xm, ym.
Тогда 1. 9 можно переписать в виде 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=w№0. тогда a1=1-w, b1=b2=1/2w и соотношения 1. 9, 1. 10, 1. 11 сведутся к
ym+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, отличных от нуля, ошибка ограничения равна
et=kh3 1. 13
Методы Рунге-Кутта третьего и четвертого порядков можно вывести совершенно аналогично тому, как это делалось при выводе методов первого и второго порядков. Мы не будем воспроизводить выкладки, а ограничимся тем, что приведем формулы, описывающие метод четвертого порядка, один из самых употребляемых методов интегрирования дифференциальных уравнений. Этот классический метод Рунге-Кутта описывается системой следующих пяти соотношений
ym+1=ym+h/6(R1+2R2+2R3+R4) 1. 14 где R1=f(xm, ym), 1. 15 R2=f(xm+h/2, ym+hR1/2), 1. 16 R3=f(xm+h/2, ym+hR2/2), 1. 17 R4=f(xm+h/2, ym+hR3/2). 1. 18
Ошибка ограничения для этого метода равна et=kh5
так что формулы 1. 14-1. 18 описывают метод четвертого порядка. Заметим, что при использовании этого метода функцию необходимо вычислять четыре раза.
3. Выбор метода реализации программы
Исходя из вышеизложенного, для решения систем дифференциальных уравнений мы выбираем наиболее точный метод решения–метод Рунге-Кутта 4 порядка, один из самых употребляемых методов интегрирования дифференциальных уравнений.
этот метод является одноступенчатым и одношаговым
требует информацию только об одной точке
имеет небольшую погрешность
значение функции рассчитывается при каждом шаге
4. Блок-схема программмы
Основная программа
Процедура INIT
Вход
f1, C[1], C[2], C[3]
f1, k1, k2, k3, k4
f1, Xn, Xk, dp, n, eps, p
выход
5. Программа
PROGRAM smith_04; USES crt; VAR i, n: integer; sum, k1, k2, k3, k4, p, dp, eps, Xn, Xk, X, dX: real; rSR, C, dC, r1, r2, r3, r4, cPR: array[1...3] of real;
f1, f2: text;
PROCEDURE Difur;
BEGIN
dC[1]: =C[3]*k2+C[2]*k4-C[1]*k1-C[1]*k3; {dcA}
dC[2]: =C[1]*k3-C[2]*k4; {dcB}
dC[3]: =C[1]*k1-C[3]*k2; {dcC}
END;
PROCEDURE RK_4;
BEGIN
Difur;
FOR i: =1 TO n DO BEGIN
r1[i]: =dC[i];
C[i]: =cPR[i]+r1[i]*(dX/2);
END;
Difur;
FOR i: =1 TO n DO BEGIN
r2[i]: =dC[i];
C[i]: =cPr[i]+r2[i]*(dX/2);
END;
Difur;
FOR i: =1 TO n DO BEGIN
r3[i]: =dC[i];
C[i]: =cPR[I]+r3[i]*dX;
END;
Difur;
FOR i: =1 TO n DO r4[i]: =dC[i];
FOR i: =1 TO n DO rSR[i]: =((r1[i]+r2[i])*(r2[i]+r3[i])*(r3[i]+r4[i]))/6; END;
PROCEDURE STROKA;
BEGIN
WRITE(f2, '|', x: 4: 1, '|', c[1]: 7: 3, '|', c[2]: 7: 3, '|', c[3]: 7: 3, '|'); WRITE(f2, sum: 3: 0, '|', dc[1]: 7: 3, '|', dc[2]: 7: 3, '|', dc[3]: 7: 3, '|'); WRITELN(f2);
END;
PROCEDURE RUN;
BEGIN
WRITE('Step 3: Calculating data and writting results to file : out. rez'); X: =Xn;
dX: =0. 05;
REPEAT
IF (ABS(x-p) Difur;
sum: =C[1]+C[2]+C[3];
STROKA;
p: =p+dp; END;
FOR i: =1 TO n DO Cpr[i]: =C[i];
RK_4;
X: =X+dX;
UNTIL(X>Xk);
WRITELN(' - done. ');
END;
PROCEDURE INIT;
BEGIN
ClrScr;
WRITELN('Smith-04: v1. 0 (c) 1998 by Mike Smith smith01@home. bar. ru '); WRITELN;
WRITELN;
WRITE('Step 1: Read data from file : in. dat'); ASSIGN(f1, 'in. dat');
RESET(f1);
READLN(f1, C[1], C[2], C[3]);
READLN(f1, k1, k2, k3, k4);
READLN(f1, Xn, Xk, dp, n, eps, p);
WRITELN(' - done. ');
ASSIGN(f2, 'out. rez');
REWRITE(f2);
WRITE('Step 2: Write header to file : out. rez'); WRITELN(f2, '=========================================================='); WRITELN(f2, '| t, c| Ca, % | Cb, %| Cc, % | SUM | dCa | dCb | dCc |'); WRITELN(f2, '=========================================================='); WRITELN(' - done. ');
END;
PROCEDURE DONE;
BEGIN
WRITELN('Step 4: Close all files and exiting.... ');
CLOSE(f1);
WRITELN(f2, '============================================================'); CLOSE(f2);
WRITELN;
END;
BEGIN
INIT;
RUN;
DONE;
END.
6. Идентификация переменных
Таблица 1
7. Результаты расчета
Таблица 2
8. Обсуждение результатов расчета.
В результате расчета кинетической схемы процесса на языке Паскаль методом Рунге-Кутты, были получены результаты зависимости изменения концентрации реагирующих веществ во времени. Исходя из полученных результатов, можно сделать вывод, что расчет произведенверно, так как, исходя из полученных значений скоростей реакций можно сделать вывод, что соблюдается баланс скоростей химической реакции.
Рассмотрим процесс подробнее. Вещество А на протяжении всего процесса расходуется на образование веществ В и С. Концентрации вещества А в начальный момент времени расходуется быстрее, чем концентрации его же в конце процесса. Это обусловлено тем, что скорость химической реакции зависит от концентрации реагирующего вещества. Производная имеет знак“минус”. Это говорит о том, что вещество расходуется. Следовательно, чем выше концентрация вещества, вступающего в процесс, тем выше скорость его реагирования с другими веществами. Вещества В и С образуются пропорционально, так как, исходя из кинетической схемы процесса и значений констант скоростей химической реакции, видно, что образование этих веществ и расходование этих веществ, одинаково. Производная имеет знак “плюс”. Это говорит о том, что вещество образуется. График. 4
Это видно также и по результатам расчета, на протяжении всего времени исследования процесса концентрации и скорости веществ В и С одинаковы. В этом можно убедиться по виду графической зависимости концентрации веществ В и С от времени.
Можно сказать, что процесс протекает в сторону увеличения концентрации веществ В и С и уменьшения концентрации веществаА. Процесс будет протекать до момента установления равновесия, но в данном случае равновесие не установлено, так как вещества продолжают расходоваться и образовываться. На протяжении всего процесса ни одно из образующихся веществ не поменяло знак производной. Это говорит о том, чтопроцесс протекает в одну сторону.
9. Инструкция к программме
Итак, программа состоит из 3 основных процедур:
Init - процедура инициалиации, включающую в себя ввод данных;
Run - процедура вычисления и обработки результатов, включает в себя вызов двух вспомогательных процедурDifur, RK-4, Stroka, первая из которых отвечает за вычисление, а последняя - за вывод результатов в файл в табличном виде;
Done - процедура подготовки к выходу из программы;
и трех вспомогательных:
Difur - процедура вычисления производных (изменение концентрации веществ за единикцу времени )
RK-4 - используя значения производных, вычисленных процедурой Difur, вычисляет последущие концентрации веществ методом Рунге-Кутта
Stroka - процедура вывода результата в файл в табличном виде
Рассмотрим все эти процедуры поподробнее:
Процедура INIT:
В данной процедуре задействованы операторы ввода/вывода Wite/Read, оператор модуля Crt - CrlScr - очистка экрана, файлового ввода/вывода - Reset/Rewrite –открытие файла для чтения и создание нового файла, соответственно. Данная процедура выполняет функцию инициализации программных данных, считывание данных из файлаin. dat, создание, открытие на запись файла out. rez и запись в него шапки таблицы результатов.
Процедура RUN:
В данной процедуре задействованы операторы цикла Repeat/Until, и For/Do c операторами условного перехода IF/Then. В зависимости от условий вызываются процедуры Difur и Strok. В теле цикла постоянно вызывается процедура RK-4 вызывающая 4 раза функциюDifur.
Процедура DONE:
В данной процедуре задействованы оператор работы с файлами Close, который закрывает файлы с исходными данными и файл с полученными в резуультате вычислений результатами.
Процедура DIFUR:
Данная процедура вычисляет производную изменения концентрации везества за единицу времени.
Процедура STROKA:
Данная процедура с помощью оператора вывода WRITEзаписывает результаты в файл, соответствующий файловой переменной F2, назначенной коммандойASSIGN в процедуре INIT
Процедура RK-4:
Данная процедура, используя вызовы процедур Difur, а также циклы операторы цикла FOR, вычисляет последуущие концентрации веществ по предидущим точкам.
Программа представляет собой 2 файла – файл с исходным текстом на языке Паскаль smith. pas и исполняемый модуль smith. exe скомпилированный компилятором TNT Pascal 3. 25 фирмы Layer`s Ins.
Исполняемый модуль программы предназначен для запуска в операционных системах: MS Dos, Windows95, Windows NT, OS/2, а также в X-windows под Linux (при наличии эмулятора )
Для нормальной работы программе необходимо 640 кb “нижней” памяти и 20 kb дискового пространства. Согласитесь –требования минимальные, учитывая то, что сама программа абсолютно не требовательна к процессору.
В процессе работы программа считывает данные из файла in. dat и записывает результаты работы в файл out. rezв табличном виде. Исходный файл программма открывает стандартными средствами ОС, не проверяя его наличие перед работой, поэтому, если данный файл не будет доступен в каталоге, в котором расположена программа, компилятор выдаст сообщение об ошибке. Если Вы после запуска программы увидели что-то типа“Runtime error 202 at 0000: 0A86”- это всего лишь значит, что программа не смогла найти файл с исходными данными в текущем каталоге. Если Вы забыли поместить его туда, скопируйте этот файл в каталог с программой и запустите исполняемый модуль еще раз. Если данный файл у Вас отсутствует, Вам прийдется сделать его самому. Для этого в любом текстовом редакторе наберите 3 выделенных строчки и сохраните созданный файл с именемin. dat
100 0 0
0. 2 0. 1 0. 2 0. 1
0 10 0. 5 3 0. 05 0
Создав файл и скопировав его к исполняемому модулю программы, запустите исполняемый модуль еще раз.
В процессе работы программа будет выдавать сообщения об успешном окончании каждого блока. Если все прошло нормально, то на экране своего компьютера Вы увидите следуще сообщения:
Step 1: Read data from file : in. dat - done. Step 2: Write header to file : out. rez - done. Step 3: Calculating data and writting results to file : out. rez - done. Step 4: Close all files and exiting....
Первый шаг (step1) сообщает, что данные из файла in. dat были успешно прочитаны
Второй – о том что программа успешно создала выходной файл out. rez и записала в него шапку таблицы с данными
В третьем сообщении сказано, что данные успешно посчитаны и записаны в выходной файлout. rez
Четвертое сообщение сообщает об окончании вычислений и завершении программы.
После того, как программа отработает, Вы сможете познакомится с результатами, которые были вычислены и помещены в файл результатовout. rez. Просмотрев его любой программой просмотра текстовых файлов или вывев его на печать, вы получите таблицу c результатами.
10. Заключение.
В результате выполнения расчета получена зависимость изменения концентрации вещества во времени. Из расчета следует, что на протяжении всего процесса вещество А расходовалось на образование В и С. Процесс не достиг конечного состояния (не достиг равновесия) Максимум концентрации вещества наблюдался при следующих значениях времени:
при начальном значении времени max соответствовал веществу А;
при значении времени, равном 10 часам, max соответствовал веществам B и С, однако, это не является максимумом концентрации веществ в процессе вообще, так как вещества B и С продолжают образовываться;
В ходе выполнения работы был произведен расчет системы дифференциальных уравнений методом Рунге-Кутты четвертого порядка, произведен расчет кинетической схемы процесса при изотермических условиях при данных значениях концентраций и констант скоростей. Расчет произведен с малой величиной погрешности.
Литература.
1. Мудров А. Е. Численные методы для ПЭВМ на языках Паскаль, Фортран и Бейсик. МП “Раско”, Томск, 1991 г.
Страницы: 1, 2
