Курсовая работа: Решение систем дифференциальных уравнений методом Рунге - Кутты 4 порядка

  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)<eps) THEN BEGIN

   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. Обсуждение результатов расчета.

В результате расчета кинетической схемы процесса на языке Паскаль методом Рунге-Кутты, были получены результаты зависимости изменения концентрации реагирующих веществ во времени. Исходя из полученных результатов, можно сделать вывод, что расчет произведен верно, так как, исходя из полученных значений скоростей реакций можно сделать вывод, что соблюдается баланс скоростей химической реакции.

Страницы: 1, 2, 3, 4