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

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

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

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

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

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

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

Таким чином знаходиться масив точок функції ym_n з різними кроками тобто 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;

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]+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);}

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("\t\t Модифiкований метод:\t Виправлений метод:");

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

{printf("\n x[%.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();}