Розв’язання звичайних диференціальних рівнянь за методом Рунге-Кутта з автоматичним вибором кроку

Размещено на http://www.allbest.ru/

Розв'язання звичайних диференціальних рівнянь за методом Рунге-Кутта з автоматичним вибором кроку

розв'язання рівняння автоматичний крок

Вступ

З огляду на те, що для методів Рунге-кутти не потрібно обчислювати додаткові початкові значення, ці методи займають особливе місце серед методів класичного типу. Нижче будуть розглянуті їх властивості, а також деякі обмеження, властиві цим методам.

Із збільшенням числа етапів для великих завдань, що вирішуються цими методами, виникли б труднощі з пам'яттю ЕОМ, крім того (і це важливіше), для великих завдань, як правило, завжди великі константи Ліпшиця. У загальному випадку це робить методи Рунге-кутти високого порядку не придатними для таких завдань. В усякому разі, інші методи зазвичай ефективніше і їм слід віддавати перевагу. Проте методи Рунге-кутти четвертого порядку є такими, що достатньо легко реалізовуються на ЕОМ, а наявність автоматичного вибору кроку дає можливість проводити обчислення з хорошою точністю. Тому їх доцільно застосовувати для досить широкої безлічі завдань.

Методи Рунге-кутти мають декілька вагомих достоїнств, що визначили їх популярність серед значного числа дослідників. Ці методи легко програмуються, володіють достатніми для широкого круга завдань властивостями точності і стійкості. Ці методи, як і всі однокрокові методи, є такими, що самостартующимі і дозволяють на будь-якому етапі обчислень легко змінювати крок інтеграції.

У роботі основна увага сконцентрована на питаннях точності і ефективності вирішення завдань того типу, для яких методи Рунге-кутти прийнятні. Програмна реалізація методів Рунге-кутти четвертого порядку з автоматичним вибором кроку представлена у вигляді програми, написаної на мові високого рівня Borland C++ 3.1. Програму можна запускати в середовищі MS-DOS або Windows® 95/98/me/2k/xp. Як вихід програма пише таблицю значень у файл на диск і малює графік на екрані ЕОМ.

1. Теоретична частина

1.1 Постановка задачі

Дано диференціальне рівняння і початкову умову, тобто поставлено завдання Коші:

Потрібно відшукати інтегральну криву, що задовольняє поставленому завданню Коші за допомогою методу Рунге-кутти четвертого порядку з автоматичним вибором кроку на відрізку . Задачу можна вирішити аналітично, знайшовши вирішення диференціального рівняння і підставивши в нього початкову умову, тим самим, відшукавши необхідну інтегральну криву. Але для нас інтерес представляє рішення даної задачі із застосуванням чисельного методу, а конкретніше - методу Рунге-кутти 4-го порядку з автоматичним вибором кроку, тобто чисельне рішення. Автоматичний вибір кроку - необхідна умова адекватної поведінки програми при функціях, що різко змінюються, задаючих інтегральну криву, що дозволяє відобразити всі моменти в поведінці інтегральної кривої і добитися високої точності.

1.2 Оцінка погрішності і збіжність методів Рунге-кутти

З часів роботи Лагранжа і особливо Коші всякий встановлений чисельно результат прийнято супроводжувати надійною оцінкою погрішності. Лагранж дав відомі оцінки погрішності многочленів Тейлора, а Коші вивів оцінки для погрішності методу ламаних Ейлера. Через декілька років після перших успіхів методів Рунге-кутти також прийшов до висновку, що для цих методів потрібні оцінки погрішностей.

2. Спеціальна частина

2.1 Постановка задачі

Програма для знаходження інтегральної кривої, задовольняючому поставленому завданню Коші написана на мові високого рівня Borland c++ 3.1. Програма складається з чотирьох функцій. За допомогою директив препроцесора #{define визначені максимальний крок і величина локальної максимальної погрішності, а також номер версії програми. Розглянемо докладніше роботу програми в комплексі.Функція title() призначена для друку на екрані назви програми.

Функція do_step() здійснює один крок Рунге-кутти і повертає набутого значення. Як вхідні параметри в неї передається поточне положення, значення шуканої функції, обчислене на попередньому кроці і величина кроку, з яким потрібно провести крок.

Функція f() задає праву частину диференціального рівняння, ліва частина диференціального рівняння рівна . Як аргументи функції передається і .

Функція main() - основна функція програми. У користувача запрошується крапка, починаючи з якої необхідно відобразити рішення задачі Коші, крапка, права межа інтеграції і значення в лівій крапці, через яке зобов'язана проходити шукана інтегральна крива. Після цього програма починає обчислювальний процес, виводячи набутих значень на екран у вигляді списку і в текстовий файл «rk4.txt» на диск. Після того, як буде досягнута права межа інтеграції, процес зупиниться і користувачеві буде запропоновано натиснути на будь-яку клавішу для того, щоб побудувати графік. Для побудови графіка програма перемикається в графічний режим. Графік масштабується з урахуванням того, щоб він завжди був видний на екрані незалежно від того, як високо або низько від осі абсцис він лежить. Крім того, програма постарається зробити його максимально наочним. Для цього будуть проведені пунктирні лінії, відповідні мінімальному і максимальному значенню інтегральної кривої, лівому і правому кінцям інтеграції, а також значенню інтегральної кривої у вказаній крапці . Для того, щоб користувач міг легко орієнтуватися на графіці, поряд з пунктирними лініями пишуться координатні значення з точністю до двох десяткових знаків. Як показали численні тести, проведені на комп'ютері на базі процесора Intel Pentium 4b з тактовою частотою 2.4 Ггц, побудова графіка відбувається значно швидше, ніж первинний розрахунок з виводом на екран і записом у файл. У цьому легко переконатися, якщо задати досить великий відрізок інтеграції, наприклад [-100,100].

Програма застосовує наступний метод Рунге-кутти четвертого порядку.

Для отримання точніших результатів і можливо навіть збільшення швидкості роботи в певних ситуаціях програма використовує автоматичне управління довжиною кроку. Проводяться одна ітерація з кроком, а потім дві ітерації підряд з кроком, починаючи з позиції x_cur. Різниця між першим і останнім згаданим тут результатами, узята по модулю, вважається локальною погрішністю. Якщо вона більша, ніж наперед задана величина, то крок ділиться навпіл і тіло циклу повторюється.

Для того, щоб програма могла «розгонитися» після зменшення кроку, передбачена умова збільшення довжини кроку. Воно полягає в тому, що якщо погрішність між другим з двох значень, обчислених з кроком і значенням, обчисленим з кроком, не перевершує, то крок збільшується удвічі і тіло циклу повторюється. Відзначимо, що величина кроку не може перевершувати значення MAXSTEP, яке визначається директивою препроцесора #define. Якщо жодне з двох описаних вище умов не виконується, це означає, що крок оптимальний і програма проводить обчислення значення функції із записом його у файл і відображенням на екрані.

Програма забезпечена механізмом захисту від збоїв - у випадку, якщо інтегрована функція терпить розрив (її не можна інтегрувати на даній ділянці), програма зупиняється і видається повідомлення про неможливість продовжувати. Працездатність цього механізму перевірена на деяких розривних функціях, таких як тангенс і ін.

Висновки

В роботі детально розглянутий метод Рунге-кутти четвертого порядку з автоматичним вибором довжини кроку, приведені необхідні теоретичні відомості, освітлені альтернативні методи і їх ефективність. Був розроблений алгоритм програмного модуля, що дозволяє автоматично міняти величину кроку інтеграції при рішенні завдання Коші залежно від необхідної точності, що є неодмінною вимогою, що пред'являється до всіх хороших сучасних програм даного класу, написаний додаток, вирішені приклади.

Перелік посилань

1. Амоносов а.А., Дубінський ю.А., Копченова н.В. «Обчислювальні методи для інженерів», М., Вища школа, 1994, 544с.

2. Хайрер Е., Нерсетт С., Ваннер Р. «Вирішення звичайних диференціальних рівнянь. Нежорсткі завдання», М., Мир, 1990, 512с.

3. Хол Д., Уатт Д. «Сучасні чисельні методи вирішення звичайних диференціальних рівнянь», М., Мир, 1979, 312с.

Додаток

Текст програми

#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();

}

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] - межі інтеграції, у(%lg)=%lg - початкова умова.\n", а, b, а, 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(" Вирішення диференціальних рівнянь методом Рунге-кутти з автоматичним вибором кроку\n");

printf(" з автоматичним вибором довжини кроку\n");

printf(" ______________________________%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);

}

Размещено на Allbest.ru