Координати точки (xk+1,
) визначають як точку перетину дотичної у=уk+f(xk,yk)(x-xk) до графіка інтегральної кривої задачі (1.1)-(1.2) в точці (xk,yk) з прямою х=хk [2].2. Розробка алгоритму розв’язання задачі
Стандартний спосіб розв’язання задачі Коші чисельними однокроковими методами – це зведення диференціальних рівнянь n-го порядку до систем з n рівнянь 1-го порядку і подальшого розв’язання цієї системи стандартними однокроковими методами:
Для рівняння
введемо заміну 
тоді для даного рівняння можна записати еквівалентну систему із двох рівнянь: 
Запишемо для кожного з цих рівнянь ітераційне рівняння:
для модифікованого методу :Ейлера:
для виправленого методу Ейлера:
Таким чином знаходиться масив точок функції ymn з різними кроками тобто n1=(b-a)/0,1=10+1 раз з кроком 0,1 і n2=(b-a)/0,05 раз з кроком 0,05. Це необхідно для оперативного визначення похибки за методом Рунге (екстраполяції Рідчардсона) [3].
Загальний вигляд похибки для цих двох методів
, де с визначається саме за методом Рунге 
, звідки с на кожному кроці обчислень знаходиться за формулою:
.Знаючи с можна знайти локальну похибку і просумувавши її по всьому діапазону інтегрування визначити загальну похибку обчислень.
Мовою програмування було обрано 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<stdio.h>
#include<conio.h>
#include<math.h>
#include<graphics.h>
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]+(h1/2)*(f(x1[i-1],ye1[i-1],zm1[i-1])+f(x1[i],ye1[i],ze1[i]));
ym1[i]=ym1[i-1]+(h1/2)*(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;
), то це точка перетину дотичної до інтегральної кривої задачі (1.1)-(1.2) в точці (хk,yk) з прямою 
Похибка на кожному кроці має порядок О(h3).
(1.11)
(1.12) 
(1.13)
у точці хk+1, яке потім уточнюють за формулою (1.13). Похибка методу на кожному кроці має порядок О(h3).
Тобто ця пряма утворює з додатним напрямком осі Ох кут 
