19335-1 (Решение систем дифференциальных уравнений методом Рунге - Кутты 4 порядка), страница 2
Описание файла
Документ из архива "Решение систем дифференциальных уравнений методом Рунге - Кутты 4 порядка", который расположен в категории "". Всё это находится в предмете "математика" из , которые можно найти в файловом архиве . Не смотря на прямую связь этого архива с , его также можно найти и в других разделах. Архив можно найти в разделе "курсовые/домашние работы", в предмете "математика" в общих файлах.
Онлайн просмотр документа "19335-1"
Текст 2 страницы из документа "19335-1"
Формулы 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 описывают метод четвертого порядка Заметим что при использовании этого метода функцию необходимо вычислять четыре раза
3. Выбор метода реализации программы
Исходя из вышеизложенного, для решения систем дифференциальных уравнений мы выбираем наиболее точный метод решения – метод Рунге-Кутта 4 порядка, один из самых употребляемых методов интегрирования дифференциальных уравнений
этот метод является одноступенчатым и одношаговым
требует информацию только об одной точке
имеет небольшую погрешность
значение функции рассчитывается при каждом шаге
4. Блок-схема программмы
О
НАЧАЛО
INIT
RUN
КОНЕЦ
сновная программа
П роцедура 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. Идентификация переменных 7. Результаты расчета Таблица 2 8 . Обсуждение результатов расчета. В результате расчета кинетической схемы процесса на языке Паскаль методом Рунге-Кутты, были получены результаты зависимости изменения концентрации реагирующих веществ во времени. Исходя из полученных результатов, можно сделать вывод, что расчет произведен верно, так как, исходя из полученных значений скоростей реакций можно сделать вывод, что соблюдается баланс скоростей химической реакции. Рассмотрим процесс подробнее. Вещество А на протяжении всего процесса расходуется на образование веществ В и С. Концентрации вещества А в начальный момент времени расходуется быстрее, чем концентрации его же в конце процесса. Это обусловлено тем, что скорость химической реакции зависит от концентрации реагирующего вещества. Производная имеет знак «минус». Это говорит о том, что вещество расходуется. Следовательно, чем выше концентрация вещества, вступающего в процесс, тем выше скорость его реагирования с другими веществами. Вещества В и С образуются пропорционально, так как, исходя из кинетической схемы процесса и значений констант скоростей химической реакции, видно, что образование этих веществ и расходование этих веществ, одинаково. Производная имеет знак «плюс» . Это говорит о том, что вещество образуется. График. 4 Это видно также и по результатам расчета, на протяжении всего времени исследования процесса концентрации и скорости веществ В и С одинаковы. В этом можно убедиться по виду графической зависимости концентрации веществ В и С от времени. Можно сказать, что процесс протекает в сторону увеличения концентрации веществ В и С и уменьшения концентрации вещества А. Процесс будет протекать до момента установления равновесия, но в данном случае равновесие не установлено, так как вещества продолжают расходоваться и образовываться. На протяжении всего процесса ни одно из образующихся веществ не поменяло знак производной. Это говорит о том, что процесс протекает в одну сторону. 9. Инструкция к программме Итак, программа состоит из 3 основных процедур: Init - процедура инициалиации, включающую в себя ввод данных; Run - процедура вычисления и обработки результатов, включает в себя вызов двух вспомогательных процедур Difur, RK-4, Stroka, первая из которых отвечает за вычисление, а последняя - за вывод результатов в файл в табличном виде; Done - процедура подготовки к выходу из программы; и трех вспомогательных: Difur - процедура вычисления производных (изменение концентрации веществ за единикцу времени ) RK-4 - используя значения производных, вычисленных процедурой Difur, вычисляет последущие концентрации веществ методом Рунге-Кутта Stroka - процедура вывода результата в файл в табличном виде Рассмотрим все эти процедуры поподробнее: Процедура INIT: PROCEDURE INIT; BEGIN ClrScr; 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; Таблица 1
В данной процедуре задействованы операторы ввода/вывода Wite/Read, оператор модуля Crt - CrlScr - очистка экрана, файлового ввода/вывода - Reset/Rewrite – открытие файла для чтения и создание нового файла, соответственно. Данная процедура выполняет функцию инициализации программных данных, считывание данных из файла in.dat, создание, открытие на запись файла out.rez и запись в него шапки таблицы результатов.
Процедура RUN:
PROCEDURE RUN;
BEGIN
X:=Xn;
dX:=0.05;
REPEAT
IF (ABS(x-p) Difur; sum:=C[1]+C[2]+C[3];