}
textcolor(0);
setbkcolor(0);
title();
printf("y'=f(x,y), y(x0)=y(a)=y0, [a,b] - отрезок интегрирования\n");
label1: printf("\na=");
scanf("%lg", &a);
printf("b=");
scanf("%lg", &b);
// Авто смена границ при необходимости
if (a > b)
{
temp = a;
a = b;
b = temp;
}
if (a == b)
{
printf("Начало отрезка интегрирования совпадает с его концом, повторите ввод!\n");
goto label1;
}
printf("y(%lg)=", a);
scanf("%lg", &y0);
title();
printf("[%lg,%lg] - границы интегрирования, y(%lg)=%lg - начальное условие.\n", a, b, a, y0);
// Инициализация
h = fabs(b - a) / 10;
if (h > 0.1) h = 0.1;
x_cur = a;
y_cur = y0;
f_max = y_cur;
f_min = y_cur;
myfile = fopen("rk4.txt", "w");
fprintf(myfile, "Program: Ilya RK4 Version %g\n", VERSION);
fprintf(myfile, "Method: Runge-Kutta\n");
fprintf(myfile, "The order of method: 4\n");
fprintf(myfile, "Automatic integration step select: Enabled\n");
fprintf(myfile, "[a,b]=[%lg,%lg], y(%lg)=%lg\n", a, b, a, y0);
while (x_cur <= b)
{
if (flag > 1) break;
big_step_res = do_step(h, x_cur, y_cur);
temp = do_step(h / 2, x_cur, y_cur);
small_step_res = do_step(h / 2, x_cur + h / 2, temp);
err = fabs(big_step_res - small_step_res);
// Уменьшение длины шага
if (err > EPSILON)
{
h = h / 2;
continue;
}
// Увеличение длины шага
big2_step_res = do_step(h, x_cur + h, big_step_res);
super_step_res = do_step(2 * h, x_cur, y_cur);
if (fabs(big2_step_res - super_step_res) < EPSILON / 2)
{
h *= 2;
continue;
}
if (h > MAXSTEP) h = MAXSTEP;
// Защита от сбоев
if (h < pow(EPSILON, 2))
{
printf("Ошибка! Возможно, функция разрывна.\nПроинтегрировать на данном интервале невозможно. Скорее всего, g(%lg)=", x_cur);
fprintf(myfile, "Ошибка! Возможно, функция разрывна.\nПроинтегрировать на данном интервале невозможно. Скорее всего, g(%lg)=", x_cur);
if (y_cur < 0)
{
printf("-oo.\n");
fprintf(myfile, "-oo.\n");
}
else
{
printf("+oo.\n");
fprintf(myfile, "+oo.\n");
}
getch();
fclose(myfile);
exit(1);
}
printf("y(%lg)=%lg, err=%lg, h=%lg\n", x_cur, y_cur, err, h);
if (y_cur < f_min) f_min = y_cur;
if (y_cur > f_max) f_max = y_cur;
fprintf(myfile, "y(%lg)=%lg, h=%lg\n", x_cur, y_cur, h);
if (x_cur + h > b) h = fabs(b - x_cur);
x_cur += h;
y_cur = big_step_res;
if (x_cur >= b) flag++;
}
fclose(myfile);
printf("\nТаблица значений записана в файл rk4.txt.\n");
printf("\nНажмите любую клавишу для построения графика...");
flag = 0;
getch();
// Построение графика
cleardevice(); clrscr();
if (fabs(a) > fabs(b)) zoom = fabs(getmaxx() / 2 / a);
else zoom = fabs(getmaxx() / 2 / b);
// Рисуем границы
for (i = 0 ; i < getmaxy() ; i += 5)
{
if (c == 8) c = 0;
else c = 8;
setcolor(c);
line(a * zoom + getmaxx() / 2, i, a * zoom + getmaxx() / 2, i + 5);
line(b * zoom + getmaxx() / 2 - 1, i, b * zoom + getmaxx() / 2 - 1, i + 5);
}
if (fabs(f_min) > fabs(f_max)) norma = fabs(f_min) * zoom;
else norma = fabs(f_max) * zoom;
// Определение коэффициента коррекции
k = (getmaxy() / 2) / norma;
// Предотвращение чрезмерного масштабирования
if (k < 0.0001) k = 0.0001;
if (k > 10000) k = 10000;
for (i = 0 ; i < getmaxx() ; i += 5)
{
if (c == 8) c = 0;
else c = 8;
setcolor(c);
line(i, -y0 * zoom * k + getmaxy() / 2, i + 5, -y0 * zoom * k + getmaxy() / 2);
line(i, -f_min * zoom * k + getmaxy() / 2, i + 5, -f_min * zoom * k + getmaxy() / 2);
line(i, -f_max * zoom * k + getmaxy() / 2, i + 5, -f_max * zoom * k + getmaxy() / 2);
}
metka = ceil((-y0 * zoom * k + getmaxy() / 2) / 16);
if (metka <= 0) metka = 1;
if (metka == 15) metka = 16;
if (metka > 25) metka = 25;
gotoxy(1, metka);
printf("Y=%.2g", y0, metka);
metka = ceil((-f_max * zoom * k + getmaxy() / 2) / 16);
if (metka <= 0) metka = 1;
if (metka == 15) metka = 16;
if (metka > 25) metka = 25;
gotoxy(1, metka);
printf("Y=%.2lg", f_max, metka);
metka = ceil((-f_min * zoom * k + getmaxy() / 2) / 16);
if (metka <= 0) metka = 1;
if (metka == 15) metka = 16;
if (metka > 25) metka = 25;
gotoxy(1, metka);
printf("Y=%.2lg", f_min, metka);
// Пишем границы, делаем отметки на осях координат
metka1 = ceil((a * zoom + getmaxx() / 2) / 8);
if (metka1 < 1) metka1 = 1;
if (metka1 > 75) metka1 = 75;
if (metka == 17) metka = 18;
gotoxy(metka1, 15);
if (a != 0) printf("%.2lg", a);
metka2 = ceil((b * zoom + getmaxx() / 2 - 1) / 8);
if (metka2 - metka1 < 7) metka2 = metka1 + 7;
if (metka2 < 1) metka2 = 1;
if (metka2 > 75) metka2 = 75;
gotoxy(metka2, 15);
printf("%.2lg", b);
gotoxy(80, 17);
printf("X");
gotoxy(42,1);
printf("Y");
gotoxy(39, 15);
printf("0");
// Рисуем систему координат
setcolor(15);
line(0, getmaxy() / 2, getmaxx(), getmaxy() / 2);
line(getmaxx() / 2, 0, getmaxx() / 2, getmaxy());
line(getmaxx() / 2, 0, getmaxx() / 2 - 5, 10);
line(getmaxx() / 2, 0, getmaxx() / 2 + 5, 10);
line(getmaxx(), getmaxy() / 2, getmaxx() - 10, getmaxy() / 2 + 5);
line(getmaxx(), getmaxy() / 2, getmaxx() - 10, getmaxy() / 2 - 5);
setcolor(10);
h = fabs(b - a) / 10;
if (h > 0.1) h = 0.1;
y_cur = y0;
x_cur = a;
f_max = y_cur;
f_min = y_cur;
x0 = zoom * a + getmaxx() / 2;
y0 = (zoom * (-y_cur)) * k + getmaxy() / 2;
while (x_cur <= b)
{
if (flag > 1) break;
big_step_res = do_step(h, x_cur, y_cur);
temp = do_step(h / 2, x_cur, y_cur);
small_step_res = do_step(h / 2, x_cur + h / 2, temp);
err = fabs(big_step_res - small_step_res);
if (err > EPSILON)
{
h = h / 2;
continue;
}
big2_step_res = do_step(h, x_cur + h, big_step_res);
super_step_res = do_step(2 * h, x_cur, y_cur);
if (fabs(big2_step_res - super_step_res) < EPSILON / 2)
{
h *= 2;
continue;
}
if (h > MAXSTEP) h = MAXSTEP;
line (x0, y0, zoom * x_cur + getmaxx() / 2, zoom * (-y_cur) * k + getmaxy() / 2);
x0 = zoom * (x_cur) + getmaxx() / 2;
y0 = (zoom * (-y_cur)) * k + getmaxy() / 2;
if (x_cur + h > b) h = fabs(b - x_cur);
x_cur += h;
y_cur = big_step_res;
if (x_cur >= b) flag++;
}
while (getch() != 0);
}
// ----------------------------------------------------------------------- //
void title(void)
{
// Печать заголовка программы
cleardevice(); clrscr();
printf(" Решение дифференциальных уравнений методом Рунге-Кутты 4-го порядка\n");
printf(" с автоматическим выбором длины шага\n");
printf(" Разработал Щербаков Илья, гр. 520212, версия %g\n", VERSION);
printf("____________________________________________________\n");
}
// ----------------------------------------------------------------------- //
double do_step(double h, double x_cur, double y_cur)
{
double k1, k2, k3, k4, delta_y_cur;
k1 = f(x_cur, y_cur);
k2 = f(x_cur + (h / 2), y_cur + (h / 2) * k1);
k3 = f(x_cur + (h / 2), y_cur + (h / 2) * k2);
k4 = f(x_cur + h, y_cur + h * k3);
delta_y_cur = (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
return(y_cur + delta_y_cur);
}
// ----------------------------------------------------------------------- //
[1] Дж. Холл, Дж. Уатт «Современные численные методы решения обыкновенных дифференциальных уравнений», М., Мир, 1979, стр. 77.
[2] «Между тем еще нет доказательства, что эти приближенные методы сходятся, или, что практически важнее, нет критерия, определяющего, сколь малым надо сделать шаги, чтобы достичь предписанной точности» – так писал Рунге в 1905 году.
[3] Хайрер Э., Нёрсетт С., Ваннер Г. «Решение обыкновенных дифференциальных уравнений. Нежесткие задачи», М., Мир, 1990, стр. 169.
[4] Амоносов А.А., Дубинский Ю.А., Копченова Н.В. «Вычислительные методы для инженеров», М., Высшая школа, 1994, стр. 445.
[5] Тестирование проводилось на компьютере на базе процессора Intel Pentium 4B. На компьютерах, оснащенных другими процессорами, время выполнения первого и второго этапов может быть другим.