8. Метод Рунге-Кутта второго порядка (Метод Эйлера-Коши)

Рассмотрим метод Рунге-Кутта второго порядка. В этом методе величины y_i+1 вычисляются по следующим формулам:

y_i+1 = y_i +  y_i

 y_i= y_i1+ y_i2

,

См. рис. 1

Тогда .

Обозначим , тогда

Геометрически это означает, что мы определяем направление интегральной кривой в исходной точке (x_i, y_i) и во вспомогательной точке , а в качестве окончательного берем среднее из этих направлений.

9. Метод Рунге-Кутта четвертого порядка

В вычислительной практике наиболее часто используется метод Рунге-Кутта четвертого порядка.

В этом методе величины y_i+1 вычисляются по следующим формулам:

y_i+1 = y_i +  y_i

 y_i=h (k₁ + 2k₂ + 2k₃ + k₄)/6, i = 0, 1, ...

k₁ = f(x_i, y_i),

k₂ = f(x_i+h/2, y_i+h k₁/2),

k₃ = f(x_i+h/2, y_i+h k₂/2),

k₄ = f(x_i+h, y_i+h k₃).

Практическая часть

2.Численная реализация решения систем дифференциальных уравнений (2) и (3)

2.1. Реализация решения в пакете MathCad методом Эйлера (3 модификация).

E0:=15

t0:=0 t1:=0.01 t2:=0.02

Графики:

1.2. Реализация решения в пакете MathCad методом Рунге-Кутта.

Графики зависимости I(t) и U(t)

1.3. Реализация решения на языке программирования высокого уровня (C++) методом Эйлера (3 модификация).

Блок-схема

Программный код

#include <iostream.h>

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

const double PI = 3.1415926;

double f, PhiA, Phi, E0, LA, L, CA, C, t0, t1, t2, h;

double R[7];

int n, i;

//Функция гармонического источника тока, учитывающая переключение ключа

double E(double t)

{

double ResultE;

if (t < t1)

{

ResultE = E0 + E0 * sin(2 * PI * f* t + Phi);

}

else

{

ResultE = 0;

}

return ResultE;

}

//Производная от тока I по времени t

double FDI(double t, double I, double U)

{

double A, B, D, ResultFDI;

A = R[2] / (R[1] + R[2]);

B = (R[1] * R[2]) / (R[1] + R[2]);

D = (R[5] + R[6]) / (R[3] + R[5] + R[6]);

ResultFDI = (1 / L) * ((A * E(t)) - (I * (R[4] + (R[3] * D) + B)) - (U * D));

return ResultFDI;

}

//Производная от напряжения U по времени t

double FDU(double t, double I, double U)

{

double A, B, ResultFDU;

A = (R[5] + R[6]) / (R[3] + R[5] + R[6]);

B = 1 / (R[3] + R[5] + R[6]);

ResultFDU = (1 / C) * ((I * A) - (U * B));

return ResultFDU;

}

void main()

{

FILE *FOUT;

double xt1, xt2, pI, I1, I2, pU, U1, U2;

FOUT = fopen("C:\\student\\Result08.txt", "a");

//Ввод параметров цепи

cout << "********** METOg 3uLEPA (3 mog.) **********\n";

cout << "BBEguTE E0 (B):\n";

cin >> E0;

cout << "BBEguTE L (MuLu rEHPu):\n";

cin >> LA;

L = LA / 1000;

cout << "BBEguTE C (MuKPO qpAPAgbI):\n";

cin >> CA;

C = CA / 1000000;

cout << "BBEguTE t0 (c):\n";

cin >> t0;

cout << "BBEguTE t1 (c);\n";

cin >> t1;

cout << "BBEguTE t2 (c):\n";

cin >> t2;

cout << "BBEguTE n:\n";

cin >> n;

h = (t2 - t0) / n;

cout << "War h = " << h << "\n";

for (i=1; i<=6; i++)

{

cout << "BBEguTE R(" << i << ") (OM):\n";

cin >> R[i];

}

cout << "BBEguTE f (Gertz):\n";

cin >> f;

cout << "BBEguTE Phi (* PI PAguAH):\n";

cin >> PhiA;

Phi = PhiA * PI;

//Определение начальных условий

cout << "BBEguTE I(0), (A):\n";

cin >> I1;

cout << "BBEguTE U(0), (B):\n";

cin >> U1;

xt1 = t0;

cout << "\n";

//Запись введенных данных в файл

fprintf(FOUT, "********** Метод Эйлера (3 мод.) **********\n");

fprintf(FOUT, "Параметры цепи:\n");

fprintf(FOUT, "E0 = %g (B)\n", E0);

fprintf(FOUT, "L = %g (мГн)\n", LA);

fprintf(FOUT, "C = %g (мкФ)\n", CA);

fprintf(FOUT, "t0 = %g (c)\nt1 = %g (c)\nt2 = %g (c)\n", t0, t1, t2);

fprintf(FOUT, "n = %d\n", n);

fprintf(FOUT, "h = %g\n", h);

for (i=1; i<=6; i++)

{

fprintf(FOUT, "R(%d) = %g (Ом)\n", i, R[i]);

}

fprintf(FOUT, "f = %g (Гц)\n", f);

fprintf(FOUT, "ф = %g*п = %g (рад.)\n", PhiA, Phi);

fprintf(FOUT, "\nТаблица значений:\n");

fprintf(FOUT, "i\tt\tI\tU\n");

fprintf(FOUT, "0\t%g\t%g\t%g\n", xt1, I1, U1);

//Цикл расчета значений функций 3 модификацией метода Эйлера

for (i=1; i<=n; i++)

{

//Вычисление значений

xt2 = xt1 + h;

pI = I1 + h * FDI(xt1, I1, U1);

pU = U1 + h * FDU(xt1, I1, U1);

I2 = I1 + (h / 2) * (FDI(xt1, I1, U1) + FDI(xt2, pI, pU));

U2 = U1 + (h / 2) * (FDU(xt1, I1, U1) + FDU(xt2, pI, pU));

//Вывод полученных значений на экран

cout << "i = " << i << "\tt = " << xt2 << "\tI = " << I2 << "\tU = " << U2 << "\n";

//Вывод полученных значений в файл

fprintf(FOUT, "%d\t%g\t%g\t%g\n", i, xt2, I2, U2);

//Смещение значений для следующей итерации

xt1 = xt2;

I1 = I2;

U1 = U2;

}

cout << "\n\n\n";

fprintf(FOUT, "\nРасчет окончен\n\n\n");

fclose(FOUT);

}

Результаты работы программы (содержание файла данных см.Приложение 1)

Сравнение результатов численной реализации решения систем дифференциальных уравнений (2) и (3):

1. сравнение результатов алгоритмов модифицированного метода Эйлера и метода Рунге-Кутта

Результаты практически идентичны друг другу.