Програма розв"язання звичайних диференціальних рівнянь однокроковими методами

Зміст

1. Короткі теоретичні відомості

2. Розробка алгоритму розв’язання задачі

3. Результати обчислень і оцінка похибки Висновки Література Додаток

1. Короткі теоретичні відомості

Часто задачі техніки і природознавства математично зводяться до відшукання розв’язку певного диференціального рівняння (або системи таких рівнянь), який задовольняє певні початкові умови (задачі Коші). Про інтегрувати таке рівняння в скінченому вигляді вдається досить рідко. при цьому дістають здебільшого такий вигляд, до якого шукана функція входить неявно, а тому користуватись ним не зручно.

На практиці застосовують здебільшого наближене інтегрування диференціальних рівнянь. Воно дає змогу знайти наближений розв’язок задачі Коші або у вигляді певного аналітичного виразу (наприклад, ряду Тейлора), або у вигляді деякої таблиці значень.

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

Задача Коші полягає в тому, щоб знайти розв’язок y (x) диференціального рівняння

(1.1)

який задовольняє початкову умову

(1.2)

Геометрично це означає, що треба знайти ту інтегральну криву y (x) рівняння (1.1), яка проходить через точку (x₀, y₀).

Задача Коші (1.1) — (1.2) має єдиний розв’язок, наприклад при виконанні умови такої теореми.

Теорема (Пікара). Якщо функція f (x, y) двох змінних х і у неперервна в замкнутому прямокутнику з центром у точці (х_0,у₀) і задовольняє в ньому умову Лівшиця по змінній у, тобто існує число K>0, яке не залежить від х і у, таке, що

(1.3)

для будь-яких точок (х₁, у₁) і (х₂, у₂), то існує єдина диференційована функція, яка є розв’язком диференціального рівняння (1.1). Цей розв’язок визначений і неперервно диференційований принаймні на відрізку [x₀-h; x₀+h], де

(1.4)

Розглянемо так звані однокрокові чисельні методи розв’язування задачі Коші (1.1)-(1.2), в яких, щоб знайти наближений розв’язок у точці х_k+1=x_k+h, досить знайти її розвязок в точці х_k.

І оскільки розв’язок задачі в точці х₀ відомий з початкових умов, то ці методи дають змогу послідовно обчислити значення розв’язку в наступних точках х₁=х₀+h, x₂=x₁+h,…

Окремим представником однокрокових чисельних методів є методи типу Ейлера. Надалі припускатимемо, що функція f (x, y) рівняння (1.1) задовольняє умови теореми Пікара.

Метод Ейлера

Нехай на відрізку [x₀, x₀+l] треба знайти чисельний розв’язок задачі Коші(1.1)-(1.2). Для цього відрізок [x₀, x₀+l] поділимо на n (для простоти) рівних частин точками

х₀, х₁, х₂,…, х_n=x₀+l, де х_k=x₀+kh (k=0,1,2,…, n), .

Величину h називають кроком чисельного інтегрування диференціального рівняння (1.1).

Розв’язати задачу (1.1)-(1.2) чисельно — це означає для заданої послідовності х₀, х₁,…, х_n=b=x₀+l незалежної змінної х та числа у₀ знайти числову послідовність у₁, у₂,…, у_n, тобто для заданої послідовності значень незалежної змінної x_k=x₀+kh (k=0, 1, …, n) побудувати таблицю наближених значень шуканого розв’язку задачі Коші.

Якщо наближений розв’язок задачі (1.1)-(1.2) в точці х_k відомий, то, проінтегрувавши рівняння (1.1) в межах від х_k до х_k+1, знайдемо його розв’язок в точці х_k+1 за формулою:

(1.5)

Саме ця формула є вихідною для побудови багатьох чисельних методів розв’язування задачі (1.1) — (1.2). Якщо інтеграл у правій частині формули (1.5) обчислити за формулою лівих прямокутників, то знайдемо

(1.6)

Відкинувши в цій рівності доданок порядку О (h²), дістанемо розрахункову формулу:

(1.7)

яку називають формулою Ейлера. у_k i y (x_k) — відповідно наближене і точне значення шуканого розв’язку задачі (1.1) і (1.2) у точці х_k. Різницю у_k-y (x_k) називають похибкою наближеного значення у_k у точці x_k.

Оскільки дотична до графіка функція у (х) в точці (x_k, y_k) матиме вигляд:

або

Звідси для ординати точки у_k+1 перетину цієї дотичної з прямою х=х_k+1 дістанем формулу (1.7), а це означає, що на кожному з відрізків [x_k, x_k+1], (k=0, 1, 2, …, n-1) інтегральна крива наближено замінюється відрізком дотичної до неї в точці (x_k, y_k). Якщо в площині О_ху позначити точки М_k(x_k;y_k), k=0, 1, 2,…, n і сполучити їх по порядку відрізками, то дістанемо ламану (її називають ламаною Ейлера), яка наближено зображує графік шуканого розв’язку задачі (1.1) — (1.2). У цьому і полягає геометричний зміст методу Ейлера (див. рис. 1)

Зазначимо, що похибка методу Ейлера на кожному кроці є величина порядку О (h²). Точність методу досить мала і переходом від точки x_k до точки x_k+1 її похибка систематично зростає.

Виправлений метод Ейлера.

Якщо інтеграл у правій частині формули (1.5) обчислити за формулою середніх прямокутників, тобто значення підінтегральної функції f (x, y (x)) обчислити в точці

то знайдемо

(1.8)

Величину невідомого значення функції у () обчислимо за формулою (1.6) з кроком. Матимемо:

Підставивши це значення у () в (1.8), дістанемо Відкинувши тут доданок пропорційний h³, матимемо Розрахункові формули вдосконаленого методу Ейлера можна записати у вигляді

Отже, в удосконаленому методі Ейлера спочатку за метод Ейлера обчислюють наближений розв’язок у задачі (1.1)-(1.2) в точці а потім наближений розв’язок у_k+1 у точці х_k+1; на кожному кроці інтегрування праву частину рівняння (1.1) обчислюють двічі (у точках (х_k, у_k) і ()).

Геометрично це означає, що на відрізку [x_k, x_k+1] графік інтегральної кривої задачі (1.1)-(1.2) замінюється відрізком прямої, яка проходить через точку (x_k, y_k) і має кутовий коефіцієнт k=. Іншими словами, ця пряма утворює з додатним напрямом осі Ох кут .

Що ж до точки (), то це точка перетину дотичної до інтегральної кривої задачі (1.1)-(1.2) в точці (х_k, y_k) з прямою Похибка на кожному кроці має порядок О (h³).

Модифікований метод Ейлера.

Якщо інтеграл в правій частині формули (1.5)обчислити за формулою трапеції, то матимемо

(1.11)

Невідоме значення у (х_k+1), що входить до правої частини цієї рівності, можна обчислити за формулою (1.7). Підставивши його в праву частину рівності (1.11), дістанемо рівність:

Звідси для удосконаленого методу Ейлера-Коші матимемо такі розрахункові формули:

(1.12)

(1.13)

Отже, і в цьому методі на кожному кроці інтегрування праву частину рівняння (1.1) обчислюють двічі: спочатку за методом Ейлера (формула (1.12)) обчислюють наближене значення шуканого розв’язку у точці х_k+1, яке потім уточнюють за формулою (1.13). Похибка методу на кожному кроці має порядок О (h³).

Така побудова наближеного розв’язку задачі (1.1)1(1.2) з геометричної точки зору означає, що на відрізку [x_k, x_k+1] графік інтегральної кривої наближають відрізком прямої, яка проходить через точку (x_k, y_k) і має кутовий коефіцієнт Тобто ця пряма утворює з додатним напрямком осі Ох кут Координати точки (x_k+1,) визначають як точку перетину дотичної у=у_k+f (x_k, y_k)(x-x_k) до графіка інтегральної кривої задачі (1.1)-(1.2) в точці (x_k, y_k) з прямою х=х_k.

2. Розробка алгоритму розв’язання задачі

Стандартний спосіб розв’язання задачі Коші чисельними однокроковими методами — це зведення диференціальних рівнянь n-го порядку до систем з n рівнянь 1-го порядку і подальшого розв’язання цієї системи стандартними однокроковими методами:

Для рівняння введемо заміну тоді для даного рівняння можна записати еквівалентну систему із двох рівнянь:

Запишемо для кожного з цих рівнянь ітераційне рівняння:

для модифікованого методу: Ейлера:

для виправленого методу Ейлера:

Таким чином знаходиться масив точок функції ym_n з різними кроками тобто n1=(b-a)/0,1=10+1 раз з кроком 0,1 і n2=(b-a)/0,05 раз з кроком 0,05. Це необхідно для оперативного визначення похибки за методом Рунге (екстраполяції Рідчардсона).

Загальний вигляд похибки для цих двох методів, де с визначається саме за методом Рунге, звідки с на кожному кроці обчислень знаходиться за формулою:

.

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

Мовою програмування було обрано Turbo C++. Вона виявилась найзручнішою із тих мов, в яких мені доводилось працювати.

Програма складається з трьох допоміжних функцій float f (x, y, z), void eylermod () i eylerisp (). eylermоd () реалізовує модифікований метод Ейлера, eylerisp () — виправлений метод, а функція f (x, y, z) повертає значення другої похідної рівняння.

Лістинг програми приведено в додатку.

3. Результати обчислень і оцінка похибки

Результатом розв’язання задачі Коші являється функція. В даному випадку отримати цю функцію в аналітичному вигляді обчислювальні однокрокові методи не дозволяють. Вони представляють функцію в табличному вигляді, тобто набір точок значень х і відповідних їм значень функції у (х). Тому для більшої наглядності було вирішено по цим точкам намалювати графіки функцій у (х) для кожного з методів окремо (дивись рисунок 4). На тому ж малюнку виведені значення похибок для кожного методу окремо. На рисунку 5 виведено значення функції у (х) в дискретному вигляді з кроком h1=0.1.

Рисунок 4.

Рисунок 5.

Висновки

В даній курсовій роботі я ознайомився з однокроковими методами розв’язання звичайних диференціальних рівнянь. Завдяки їй я остаточно розібрався застосовуванням цих методів до розв’язання диференціальних рівнянь вищих порядків на прикладі рівняння другого порядку.

Література

1. Мудров А. Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. — Томск: МП «Раско», 1991. — 272 с.

2. Бортків А.Б., Гринчишин Я. Т. Turbo Pascal: Алгоритми і програми: чисельні методи в фізиці і математиці. Навчальний посібник. — К.: Вища школа, 1992. — 247 с.

3. Квєтний Р. Н. Методи комп’ютерних обчислень. Навчальний посібник. — Вінниця: ВДТУ, 2001 — 148 с.

Додаток

Лістинг програми

#include

#include

#include

#include

float f (float x, float y, float z)

{return 0.7*z+x*y+0.7*x;}

float h1=0.1;

float h2=0.05;

float a=0;

float b=1;

float x2[21], ye2[21], ym1[11], zm2[21], ym2[21], ye1[11];

float ze1[11], zm1[11], ze2[21], x1[11], yi1[11], yi2[21];

float zi1[11], zi2[21];

int n1=(b-a)/h1;

int n2=(b-a)/h2;

void eylermod ()

{// printf («[0] %5.2f %5.2f %5.2f», x2[0], y2[0], z2[0]);

// moveto ((x2[0])*100,480-((ym2[0])*100));

for (int i=1;i<=n2+1;i++)

{x2[i]=x2[i-1]+h2;

ze2[i]=ze2[i-1]+h2*f (x2[i-1], ye2[i-1], ze2[i-1]);

ye2[i]=ye2[i-1]+h2*ze2[i-1];

zm2[i]=zm2[i-1]+(h2/2)*(f (x2[i-1], ye2[i-1], zm2[i-1])+f (x2[i], ye2[i], ze2[i]));

ym2[i]=ym2[i-1]+(h2/2)*(ze2[i]+zm2[i-1]);

// printf («n[%d] %5.2f %5.2f %5.2f», i, x2[i], ye2[i], ym2[i]);

// setcolor (YELLOW);

// lineto ((x2[i])*100,480-((ym2[i])*100));}

moveto ((x1[0])*250+20,480-((ym1[0])*100)-30);

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

{x1[i]=x1[i-1]+h1;

ze1[i]=ze1[i-1]+h1*f (x1[i-1], ye1[i-1], ze1[i-1]);

ye1[i]=ye1[i-1]+h1*ze1[i-1];

zm1[i]=zm1[i-1]+(h½)*(f (x1[i-1], ye1[i-1], zm1[i-1])+f (x1[i], ye1[i], ze1[i]));

ym1[i]=ym1[i-1]+(h½)*(ze1[i]+zm1[i-1]);

// printf («n[%d] %5.2f %5.2f %5.2f», i, x1[i], ye1[i], ym1[i]);

setcolor (12);

lineto ((x1[i])*250+20,480-((ym1[i])*100)-30);}

float c;

float s=0;

for (i=0;i<=n1+1;i++)

{c=(ym2[i*2]-ym1[i])/(h1*h1*h1-h2*h2*h2);

s+=c*h1*h1*h1;}

char *ch;

sprintf (ch," %f", fabs (s));

setcolor (15);

settextstyle (0,0,1);

outtextxy (5,108," Похибка:");

settextstyle (2,0,5);

outtextxy (70,102,ch);}

void eylerisp ()

{// printf («[0] %5.2f %5.2f %5.2f», x2[0], y2[0], z2[0]);

// moveto ((x2[0])*100,480-((ym2[0])*100));

for (int i=1;i<=n2+1;i++)

{x2[i]=x2[i-1]+h2/2;

ze2[i]=ze2[i-1]+h2/2*f (x2[i-1], ye2[i-1], ze2[i-1]);

ye2[i]=ye2[i-1]+h2/2*ze2[i];

zi2[i]=zi2[i-1]+h2*f (x2[i], ye2[i], ze2[i]);

yi2[i]=yi2[i-1]+h2*zi2[i];

x2[i]+=h2/2;

// printf («n[%d] %5.2f %5.2f %5.2f», i, x2[i], ye2[i], ym2[i]);

// setcolor (YELLOW);

// lineto ((x2[i])*100,480-((ym2[i])*100));}

moveto ((x1[0])*250+350,480-((yi1[0])*100)-30);

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

{x1[i]=x1[i-1]+h½;

ze1[i]=ze1[i-1]+h½*f (x1[i-1], ye1[i-1], ze1[i-1]);

ye1[i]=ye1[i-1]+h½*ze1[i];

zi1[i]=zi1[i-1]+h1*f (x1[i], ye1[i], ze1[i]);

yi1[i]=yi1[i-1]+h1*zi1[i];

x1[i]+=h½;

// printf («n[%d] %5.2f %5.2f %5.2f», i, x1[i], ye1[i], ym1[i]);

setcolor (12);

lineto ((x1[i])*250+350,480-((yi1[i])*100)-30);}

float c;

float s=0;

for (i=0;i<=n1+1;i++)

{c=(yi2[i*2]-yi1[i])/(h1*h1*h1-h2*h2*h2);

s+=c*h1*h1*h1;}

char *ch;

sprintf (ch," %f", fabs (s));

setcolor (15);

settextstyle (0,0,1);

outtextxy (335,108," Похибка:");

settextstyle (2,0,5);

outtextxy (405,102,ch);}

void main ()

{float c=0,s=0;

int gdriver = DETECT, gmode, errorcode;

initgraph (&gdriver, &gmode, «»);

cleardevice ();

x2[0]=x1[0]=a;

ye2[0]=ye1[0]=1;

ym2[0]=ym1[0]=1;

ze2[0]=ze1[0]=1;

zm2[0]=zm1[0]=1;

yi2[0]=yi1[0]=1;

zi2[0]=zi1[0]=1;

char v=50;

while (v≠27)

{//setgraphmode (getgraphmode ());

setbkcolor (16);

outtextxy (190,0," Курсова робота з дисциплiни");

setcolor (10);

outtextxy (205,10," <<�Обчислювальнi методи>>");

setcolor (12);

outtextxy (95,20," на тему: <<�Дослiдження однокрокових методiв розв’язання");

outtextxy (165,30," звичайних диференцiальних рiвнянь>>");

setcolor (14);

outtextxy (25,90," Модифiкованний метод Ейлера");

outtextxy (355,90," Виправлений метод Ейлера");

setcolor (15);

outtextxy (455,50," Виконав ст. гр. 2АВ-01″);

outtextxy (455,60," Сторожук Костянтин");

settextstyle (8,0,1);

outtextxy (45,45," y''=0.7y'+xy+0.7x");

settextstyle (0,0,1);

setcolor (7);

line (20,475,20,120); //левая ось у

line (0,450,300,450); //левая ось х

line (350,475,350,120);//правая ось у

line (330,450,630,450);//правая ось х

line (20,120,18,130);

line (20,120,22,130); //стрелки оу

line (18,130,22,130);

line (300,450,290,448);

line (300,450,290,452); //срелки ох

line (290,448,290,452);

line (350,120,348,130);

line (350,120,352,130); //стрелки оу

line (348,130,352,130);

line (630,450,620,448);

line (630,450,620,452); //срелки ох

line (620,448,620,452);

char t[5];

char m[5];

settextstyle (2,0,5);

outtextxy (285,430," x");

outtextxy (28,122," y (x)");

outtextxy (615,430," x");

outtextxy (358,122," y (x)");

for (float i=0;i<11;i++)

{line (20+i*25,447,20+i*25,453);

if (i<10)line (18,450-(i*50)/1.5,22,450-(i*50)/1.5);

sprintf (t," %.1f", i/10);

if (int (i)%2==0) outtextxy (i*25+12,460,t);

sprintf (m," %.0f", i+1);

if (i<3)outtextxy (8,342-i*100,m);}

for (i=0;i<11;i++)

{line (350+i*25,447,350+i*25,453);

if (i<10)line (348,450-(i*50)/1.5,352,450-(i*50)/1.5);

sprintf (t," %.1f", i/10);

if (int (i)%2==0) outtextxy (i*25+342,460,t);

sprintf (m," %.0f", i+1);

if (i<3)outtextxy (338,342-i*100,m);}

settextstyle (0,0,1);

eylermod ();

eylerisp ();

v=getch ();

if (v==27)break;

restorecrtmode ();

setgraphmode (getgraphmode ());

printf («tt Модифiкований метод: t Виправлений метод:»);

for (i=0;i

{printf («n x[%.f]=%.1ftty (x)=%f tt y (x)=%f», i, x1[i], ym1[i], yi1[i]);}

settextstyle (0,0,1);

v=getch ();

cleardevice ();}

closegraph ();}