Зміст
1. Короткі теоретичні відомості
2. Розробка алгоритму розв’язання задачі
3.Результати обчислень і оцінка похибки
Висновки
Література
Додаток
1.Короткі теоретичні відомості
Часто задачітехніки і природознавства математично зводяться до відшукання розв’язку певногодиференціального рівняння (або системи таких рівнянь), який задовольняє певніпочаткові умови (задачі Коші). Про інтегрувати таке рівняння в скінченомувигляді вдається досить рідко. при цьому дістають здебільшого такий вигляд, доякого шукана функція входить неявно, а тому користуватись ним не зручно.
На практицізастосовують здебільшого наближене інтегрування диференціальних рівнянь. Вонодає змогу знайти наближений розв’язок задачі Коші або у вигляді певногоаналітичного виразу (наприклад, ряду Тейлора), або у вигляді деякої таблицізначень.
Розглянемо окреміметоди чисельного розв’язування задачі Коші для звичайного диференціальногорівняння першого порядку, розв’язаного відносно похідної. Наближений розв’язокзадачі Коші записують у вигляді певної таблиці значень.
Задача Кошіполягає в тому, щоб знайти розв’язок y(x) диференціального рівняння
/>, (1.1)
який задовольняєпочаткову умову
/>(1.2)
Геометрично цеозначає, що треба знайти ту інтегральну криву y(x) рівняння (1.1), якапроходить через точку (x0,y0).
Задача Коші (1.1)– (1.2) має єдиний розв’язок, наприклад при виконанні умови такої теореми.
Теорема (Пікара). Якщо функція f(x,y)двох змінних х і у неперервна в замкнутому прямокутнику
/>
з центром у точці(х0, у0) і задовольняє в ньому умову Лівшиця по змінній у,тобто існує число K>0, яке не залежить від х і у, таке, що
/>(1.3)
для будь-якихточок (х1, у1)/> і(х2, у2) />, тоіснує єдина диференційована функція />, яка єрозв’язком диференціального рівняння (1.1). Цей розв’язок визначений інеперервно диференційований принаймні на відрізку [x0-h; x0+h],де
/>(1.4)
Розглянемо такзвані однокрокові чисельні методи розв’язування задачі Коші (1.1)-(1.2), вяких, щоб знайти наближений розв’язок у точці хk+1=xk+h,досить знайти її розвязок в точці хk.
І оскількирозв’язок задачі в точці х0відомий з початкових умов, то ціметоди дають змогу послідовно обчислити значення розв’язку в наступних точках х1=х0+h,x2=x1+h,…
Окремимпредставником однокрокових чисельних методів є методи типу Ейлера. Надаліприпускатимемо, що функція f(x,y) рівняння (1.1) задовольняє умовитеореми Пікара [1].
Метод Ейлера
Нехай на відрізку[x0,x0+l] треба знайти чисельний розв’язок задачіКоші(1.1)-(1.2). Для цього відрізок [x0,x0+l]поділимо на n (для простоти) рівних частин точками
х0,х1, х2,..., хn=x0+l, де хk=x0+kh(k=0,1,2,...,n), />.
Величину hназивають кроком чисельного інтегрування диференціального рівняння (1.1).
Розв’язати задачу(1.1)-(1.2) чисельно – це означає для заданої послідовності х0, х1,…,хn=b=x0+l незалежної змінної х та числа у0знайти числову послідовність у1, у2,…, уn,тобто для заданої послідовності значень незалежної змінної xk=x0+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)
Відкинувши в ційрівності доданок порядку О(h2), дістанемо розрахункову формулу:
/> /> /> (1.7)
яку називають формулоюЕйлера. уk i y(xk) – відповідно наближене і точне значенняшуканого розв’язку задачі (1.1) і (1.2) у точці хk. Різницю уk-y(xk)називають похибкою наближеного значення уk у точці xk.
Оскільки дотичнадо графіка функція у(х) в точці (xk,yk) матиме вигляд:
/> або />
Звідси дляординати точки уk+1 перетину цієї дотичної з прямою х=хk+1дістанем формулу (1.7), а це означає, що на кожному з відрізків [xk,xk+1],(k=0, 1, 2, ..., n-1 ) інтегральна крива наближено замінюється відрізкомдотичної до неї в точці (xk,yk). Якщо в площині Охупозначити точки Мk(xk;yk), k=0, 1, 2,...,n ісполучити їх по порядку відрізками, то дістанемо ламану (її називають ламаноюЕйлера), яка наближено зображує графік шуканого розв’язку задачі (1.1) – (1.2).У цьому і полягає геометричний зміст методу Ейлера (див. рис. 1)
/>/>/>/>/>/>
/>
/>/>
/>
х4
х1
х
х0
Зазначимо, щопохибка методу Ейлера на кожному кроці є величина порядку О(h2).Точність методу досить мала і переходом від точки xk до точки xk+1її похибка систематично зростає.
Виправленийметод Ейлера.
Якщо інтеграл управій частині формули (1.5) обчислити за формулою середніх прямокутників,тобто значення підінтегральної функції f(x,y(x)) обчислити в точці
/>, то знайдемо
/> (1.8)
Величинуневідомого значення функції у(/>)обчислимо за формулою (1.6) з кроком />.Матимемо:
/>
Підставивши цезначення у(/>) в (1.8), дістанемо
/>
Відкинувши тутдоданок пропорційний h3, матимемо
/>/>
Розрахунковіформули вдосконаленого методу Ейлера можна записати у вигляді
/>
Отже, вудосконаленому методі Ейлера спочатку за метод Ейлера обчислюють наближенийрозв’язок у/> задачі (1.1)-(1.2) в точці/> а потім наближенийрозв’язок уk+1 у точці хk+1; на кожному кроціінтегрування праву частину рівняння (1.1) обчислюють двічі (у точках (хk, уk)і (/>)).
Геометрично цеозначає, що на відрізку [xk,xk+1] графік інтегральноїкривої задачі (1.1)-(1.2) замінюється відрізком прямої, яка проходить черезточку (xk,yk) і має кутовий коефіцієнт k=/>. Іншими словами, ця прямаутворює з додатним напрямом осі Ох кут />. /> /> /> /> /> /> /> />
Рисунок 2.
xk
xk+/>
xk+1
x
Y
Yk+/>
Yk
y=f(x)
l2
l1
М
Yk+1
Що ж до точки (/>), то це точка перетинудотичної до інтегральної кривої задачі (1.1)-(1.2) в точці (хk,yk)з прямою /> Похибка на кожному кроцімає порядок О(h3).
Модифікованийметод Ейлера.
Якщо інтеграл вправій частині формули (1.5)обчислити за формулою трапеції, то матимемо
/> (1.11)
Невідоме значенняу(хk+1), що входить до правої частини цієї рівності, можна обчислитиза формулою (1.7). Підставивши його в праву частину рівності (1.11), дістанеморівність:
/>
Звідси дляудосконаленого методу Ейлера-Коші матимемо такі розрахункові формули:
/> (1.12)
/> (1.13)
Отже, і в цьомуметоді на кожному кроці інтегрування праву частину рівняння (1.1) обчислюютьдвічі: спочатку за методом Ейлера (формула (1.12)) обчислюють наближенезначення шуканого розв’язку /> у точціхk+1, яке потім уточнюють за формулою (1.13). Похибка методу накожному кроці має порядок О(h3).
Така побудованаближеного розв’язку задачі (1.1)1(1.2) з геометричної точки зору означає, щона відрізку [xk,xk+1] графік інтегральної кривоїнаближають відрізком прямої, яка проходить через точку (xk,yk)і має кутовий коефіцієнт /> Тобтоця пряма утворює з додатним напрямком осі Ох кут
/>
/>
M
l1
l
l
l2
/>k+1
yk+1
yk
хk
хk+1
х
Координати точки(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
#include
#include
#include
floatf(float x,float y,float z)
{return0.7*z+x*y+0.7*x;}
floath1=0.1;
floath2=0.05;
floata=0;
floatb=1;
floatx2[21],ye2[21],ym1[11],zm2[21],ym2[21],ye1[11];
floatze1[11],zm1[11],ze2[21],x1[11],yi1[11],yi2[21];
floatzi1[11],zi2[21];
intn1=(b-a)/h1;
intn2=(b-a)/h2;
voideylermod()
{//printf("[0] %5.2f %5.2f %5.2f",x2[0],y2[0],z2[0]);
//moveto((x2[0])*100,480-((ym2[0])*100));
for(inti=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
{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);}
floatc;
floats=0;
for(i=0;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);}
voideylerisp()
{//printf("[0] %5.2f %5.2f %5.2f",x2[0],y2[0],z2[0]);
//moveto((x2[0])*100,480-((ym2[0])*100));
for(inti=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
{x1[i]=x1[i-1]+h1/2;
ze1[i]=ze1[i-1]+h1/2*f(x1[i-1],ye1[i-1],ze1[i-1]);
ye1[i]=ye1[i-1]+h1/2*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]+=h1/2;
//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);}
floatc;
floats=0;
for(i=0;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);}
voidmain()
{floatc=0,s=0;
intgdriver = 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;
charv=50;
while(v!=27)
{//setgraphmode(getgraphmode());
setbkcolor(16);
outtextxy(190,0,«Курсоваробота з дисциплiни»);
setcolor(10);
outtextxy(205,10,">");
setcolor(12);
outtextxy(95,20,«натему:
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);
chart[5];
charm[5];
settextstyle(2,0,5);
outtextxy(285,430,«x»);
outtextxy(28,122,«y(x)»);
outtextxy(615,430,«x»);
outtextxy(358,122,«y(x)»);
for(floati=0;i
{line(20+i*25,447,20+i*25,453);
if(i
sprintf(t,"%.1f",i/10);
if(int(i)%2==0)outtextxy(i*25+12,460,t);
sprintf(m,"%.0f",i+1);
if(i
for(i=0;i
{line(350+i*25,447,350+i*25,453);
if(i
sprintf(t,"%.1f",i/10);
if(int(i)%2==0)outtextxy(i*25+342,460,t);
sprintf(m,"%.0f",i+1);
if(i
settextstyle(0,0,1);
eylermod();
eylerisp();
v=getch();
if(v==27)break;
restorecrtmode();
setgraphmode(getgraphmode());
printf("\t\tМодифiкований метод:\t Виправлений метод:");
for(i=0;i
{printf("\nx[%.f]=%.1f\t\ty(x)=%f \t\t y(x)=%f",i,x1[i],ym1[i],yi1[i]);}
settextstyle(0,0,1);
v=getch();
cleardevice();}
closegraph();}