Метод Рунге-Кутты четвертого порядка с автоматическим выбором шага интегрирования решения задачи Коши (стр. 6 из 7)
Программа применяет следующий метод Рунге-Кутты четвертого порядка.
Для получения более точных результатов и возможно даже увеличения скорости работы в определенных ситуациях программа использует автоматическое управление длиной шага. Производятся одна итерация с шагом
, а затем две итерации подряд с шагом 
, начиная с позиции x_cur. Разность между первым и последним упомянутым здесь результатами, взятая по модулю, считается локальной погрешностью. Если она больше, чем наперед заданная величина 
, то шаг делится пополам и тело цикла повторяется.Для того, чтобы программа могла «разгоняться» после уменьшения шага, предусмотрено условие увеличения длины шага. Оно состоит в том, что если погрешность между вторым из двух значений, вычисленных с шагом
и значением, вычисленным с шагом 
, не превосходит 
, то шаг увеличивается вдвое и тело цикла повторяется. Отметим, что величина шага не может превосходить значения MAXSTEP, которое определяется директивой препроцессора #define. Если ни одно из двух описанных выше условий не выполняется, это означает, что шаг оптимален и программа производит вычисление значения функции с записью его в файл и отображением на экране.Программа снабжена механизмом защиты от сбоев – в случае, если интегрируемая функция терпит разрыв (ее нельзя интегрировать на данном участке), программа останавливается и выдается сообщение о невозможности продолжать. Работоспособность этого механизма проверена на некоторых разрывных функциях, таких как тангенс и др.
Заключение
В работе детально рассмотрен метод Рунге-Кутты четвертого порядка с автоматическим выбором длины шага, приведены необходимые теоретические сведения, освещены альтернативные методы и их эффективность.
Был разработан алгоритм программного модуля, позволяющий автоматически менять величину шага интегрирования при решении задачи Коши в зависимости от требуемой точности, что является непременным требованием, предъявляемым ко всем хорошим современным программам данного класса, написано приложение, решены примеры.
Список использованных источников
[1]. Амоносов А.А., Дубинский Ю.А., Копченова Н.В. «Вычислительные методы для инженеров», М., Высшая школа, 1994, 544с.
[2]. Хайрер Э., Нёрсетт С., Ваннер Г. «Решение обыкновенных дифференциальных уравнений. Нежесткие задачи», М., Мир, 1990, 512с.
[3]. Холл Д., Уатт Д. «Современные численные методы решения обыкновенных дифференциальных уравнений», М., Мир, 1979, 312с.
Приложение А. Графики функций
В данном приложении рассмотрены три дифференциальных уравнения первого порядка. К каждому уравнению прилагается по два графика – первый из них построен созданным приложением, а второй создан в пакете Maple 9.01.
Интегральная кривая, построенная приложением «Ilya RK-4 версия 1.43»
Интегральная кривая, построенная математическим пакетом Waterloo Maple 9.01
Интегральная кривая, построенная приложением «Ilya RK-4 версия 1.43»
Интегральная кривая, построенная математическим пакетом Waterloo Maple 9.01
Интегральная кривая, построенная приложением «Ilya RK-4 версия 1.43»
Интегральная кривая, построенная математическим пакетом Waterloo Maple 9.01
Интегральная кривая, построенная приложением «Ilya RK-4 версия 1.43»
Интегральная кривая, построенная математическим пакетом Waterloo Maple 9.01
Приложение Б.
Пример таблицы значений функции y(x)
y(-7)=100, h=0.4
y(-6.6)=100.045, h=0.4
y(-6.2)=100.112, h=0.4
y(-5.8)=100.212, h=0.4
y(-5.4)=100.361, h=0.4
y(-5)=100.585, h=0.4
y(-4.6)=100.919, h=0.4
y(-4.2)=101.419, h=0.4
y(-3.8)=102.17, h=0.4
y(-3.4)=103.301, h=0.2
y(-3.2)=104.067, h=0.2
y(-3)=105.011, h=0.2
y(-2.8)=106.175, h=0.2
y(-2.6)=107.615, h=0.2
y(-2.4)=109.4, h=0.2
y(-2.2)=111.62, h=0.2
y(-2)=114.392, h=0.2
y(-1.8)=117.873, h=0.2
y(-1.6)=122.267, h=0.2
y(-1.4)=127.857, h=0.2
y(-1.2)=135.033, h=0.2
y(-1)=144.346, h=0.2
y(-0.8)=156.596, h=0.2
y(-0.6)=172.977, h=0.1
y(-0.5)=183.256, h=0.1
y(-0.4)=195.327, h=0.1
y(-0.3)=209.595, h=0.1
y(-0.2)=226.578, h=0.1
y(-0.1)=246.953, h=0.05
y(-0.05)=258.68, h=0.05
y(2.96985e-15)=271.608, h=0.05
y(0.05)=285.897, h=0.05
y(0.1)=301.73, h=0.05
y(0.15)=319.32, h=0.05
y(0.2)=338.919, h=0.05
y(0.25)=360.821, h=0.05
y(0.3)=385.374, h=0.05
y(0.35)=412.989, h=0.05
y(0.4)=444.156, h=0.05
y(0.45)=479.459, h=0.025
y(0.475)=498.877, h=0.025
y(0.5)=519.603, h=0.025
y(0.525)=541.747, h=0.025
y(0.55)=565.433, h=0.025
y(0.575)=590.794, h=0.025
y(0.6)=617.978, h=0.025
y(0.625)=647.149, h=0.025
y(0.65)=678.489, h=0.025
y(0.675)=712.199, h=0.025
y(0.7)=748.501, h=0.025
y(0.725)=787.645, h=0.025
y(0.75)=829.906, h=0.025
y(0.775)=875.592, h=0.025
y(0.8)=925.047, h=0.025
y(0.825)=978.656, h=0.025
y(0.85)=1036.85, h=0.025
y(0.875)=1100.11, h=0.025
y(0.9)=1168.98, h=0.025
y(0.925)=1244.07, h=0.0125
y(0.9375)=1284.16, h=0.0125
y(0.95)=1326.08, h=0.0125
y(0.9625)=1369.91, h=0.0125
y(0.975)=1415.77, h=0.0125
y(0.9875)=1463.78, h=0.0125
y(1)=1514.04, h=0.0125
y(1.0125)=1566.7, h=0.0125
y(1.025)=1621.89, h=0.0125
y(1.0375)=1679.75, h=0.0125
y(1.05)=1740.44, h=0.0125
y(1.0625)=1804.13, h=0.0125
y(1.075)=1871, h=0.0125
y(1.0875)=1941.23, h=0.0125
y(1.1)=2015.04, h=0.0125
y(1.1125)=2092.63, h=0.0125
y(1.125)=2174.24, h=0.0125
y(1.1375)=2260.12, h=0.0125
y(1.15)=2350.54, h=0.0125
y(1.1625)=2445.79, h=0.0125
y(1.175)=2546.16, h=0.0125
y(1.1875)=2652, h=0.0125
y(1.2)=2763.65, h=0.0125
y(1.2125)=2881.49, h=0.0125
y(1.225)=3005.94, h=0.0125
y(1.2375)=3137.43, h=0.0125
y(1.25)=3276.43, h=0.0125
y(1.2625)=3423.46, h=0.0125
y(1.275)=3579.07, h=0.0125
y(1.2875)=3743.83, h=0.0125
y(1.3)=3918.41, h=0.0125
y(1.3125)=4103.47, h=0.0125
y(1.325)=4299.77, h=0.0125
y(1.3375)=4508.11, h=0.0125
y(1.35)=4729.35, h=0.00625
y(1.35625)=4845.11, h=0.00625
y(1.3625)=4964.45, h=0.00625
y(1.36875)=5087.5, h=0.00625
y(1.375)=5214.41, h=0.00625
y(1.38125)=5345.3, h=0.00625
y(1.3875)=5480.34, h=0.00625
y(1.39375)=5619.66, h=0.00625
y(1.4)=5763.44, h=0.00625
y(1.40625)=5911.83, h=0.00625
y(1.4125)=6065, h=0.00625
y(1.41875)=6223.14, h=0.00625
y(1.425)=6386.44, h=0.00625
y(1.43125)=6555.09, h=0.00625
y(1.4375)=6729.28, h=0.00625
y(1.44375)=6909.25, h=0.00625
y(1.45)=7095.2, h=0.00625
y(1.45625)=7287.36, h=0.00625
y(1.4625)=7485.99, h=0.00625
y(1.46875)=7691.33, h=0.00625
y(1.475)=7903.64, h=0.00625
y(1.48125)=8123.19, h=0.00625
y(1.4875)=8350.28, h=0.00625
y(1.49375)=8585.2, h=0.00625
y(1.5)=8828.27, h=0.00625
Приложение В.
Листинг программы «Ilya RK-4 версия 1.43»
// ----------------------------------------------------------------------- //
#include <dos.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <graphics.h>
#include <stdlib.h>
#define EPSILON 0.00001
#define MAXSTEP 1
#define VERSION 1.43
// ----------------------------------------------------------------------- //
double f(double x, double y);
double do_step(double h, double x_cur, double y_cur);
void title(void);
void main(void);
// ----------------------------------------------------------------------- //
double f(double x, double y)
{
// Правая часть ДУ f(x,y)
return (pow(2.718,x)*y);
}
// ----------------------------------------------------------------------- //
void main(void)
{
int i; // Вспомогательный счетчик
int metka; // Метка на осях
int flag = 0; // Флаг правого конца интегрирования
int metka1, metka2; // Переменные меток на осях координат
double err = 0; // Погрешность
double x0, y0; // Координаты точки начального условия
double big2_step_res, super_step_res; // Результаты длинных шагов
double k = 1; // Коэффициент коррекции
double zoom = 1; // Масштаб на графике
double big_step_res, small_step_res; // Результаты шагов интегрирования
double a, b; // Границы
double temp; // Переменная для служебных нужд
double x_cur, y_cur; // Переменные метода РК
double h; // Шаг интегрирования
double f_max = 0, f_min = 0; // Максимальное и минимальное значение кривой
double norma = 0; // Норма (для корректного масштабирования графика)int c = 8; // Переменная цвета разделительных линий
FILE *myfile; // Указатель на текстовый файл с таблицей значений
// Инициализируем графический адаптер
int gdriver = DETECT, gmode, errorcode;
initgraph(&gdriver, &gmode, "");
errorcode = graphresult();
if (errorcode != grOk)
{
printf("Ошибка инициализации графики: %s\n", grapherrormsg(errorcode));
getch();