Лабораторна робота №22
Розв’язування звичайних диференційних рівнянь методами Ейлера та Рунге-Кутта
Завдання:
Розв’язати задачу Коші (звичайне диференційне рівняння):
, початкові умови: .
1. Знайдемо аналітичний розв’язок даної задачі Коші. Аналітичний розв’язок дифрівняння в даному випадку можна знайти, використавши метод розділення змінних. Тоді отримаємо:

Визначимо із початкової умови значення сталої :

Отже, розв’язок звичайного диференційного рівняння має вигляд:

Формула методу Ейлера для числового розв’язування
звичайного диференційного рівняння

Блок-схема числового розв’язування

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

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


#include <graphics.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#define C_1 1.2*exp(cos(2))
#define min(x,y) x<y ? x : y
#define max(x,y) x>y ? x : y
float f(float x, float y)
{return y*sin(x+1);}
main()
{int gdriver = DETECT, gmode, errorcode;
/* initialize graphics mode */
initgraph(&gdriver, &gmode,"d:\\tcpp1_01\\tc\\bgi");
float x_min=1, x_max=24, y_min=1.2, y_max=1.2, sc_x, sc_y;
float x,y,y1;
float ye_0=1.2, ye_1, h;
float k1, k2, k3, k4, yrk_0=1.2, yrk_1;
float z_min, z_max;
sc_x=(x_max-x_min)/getmaxx();
for(x=x_min; x<x_max; x+=sc_x)
{y=C_1*exp(-cos(x+1));h=sc_x;
//printf("\n x=%4.2f y=%8.6f ye_0=%6.4f yrk_0=%8.6f",x,y,ye_0,yrk_0);
ye_1=ye_0+h*f(x,ye_0);
ye_0=ye_1;
k1=h*f(x,yrk_0);
k2=h*f(x+h/2,yrk_0+k1/2);
k3=h*f(x+h/2,yrk_0+k2/2);
k4=h*f(x+h,yrk_0+k3);
yrk_1=yrk_0+1.0/6*(k1+2*k2+2*k3+k4);
yrk_0=yrk_1;
z_min=min(y, min(yrk_0, ye_0));
if(y_min>z_min) y_min=z_min;
z_max=max(y, max(yrk_0, ye_0));
if(y_max<z_max) y_max=z_max;
}
//printf("\n y_min=%6.4f y_max=%6.4f",y_min,y_max);
getchar();

int x_axe, y_axe;
int xx,yy,yy1,yy2;

sc_y=(y_max-y_min)/getmaxy();
//printf("\n sc_x=%f sc_y=%f",sc_x,sc_y);
//getchar();
if(x_min*x_max<=0) x_axe=ceil(fabs(x_min)/sc_x);
if(y_min*y_max<=0) y_axe=ceil(fabs(y_max)/sc_y);
//printf("\n x_axe=%d",x_axe),printf("\n y_axe=%d",y_axe);
line(x_axe,0,x_axe,getmaxy());
line(0,y_axe,getmaxx(),y_axe);
getchar();
ye_0=1.2, yrk_0=1.2;
for(x=x_min;x<=x_max;x+=sc_x)
{y=C_1*exp(-cos(x+1)); h=sc_x;
ye_1=ye_0+h*f(x,ye_0);
ye_0=ye_1;
yy1=ye_0;
k1=h*f(x,yrk_0);
k2=h*f(x+h/2,yrk_0+k1/2);
k3=h*f(x+h/2,yrk_0+k2/2);
k4=h*f(x+h,yrk_0+k3);
yrk_1=yrk_0+1.0/6*(k1+2*k2+2*k3+k4);
yrk_0=yrk_1;
yy2=yrk_0;
//y1=2*sin(x)*sin(x);
//a*ymin+b=getmaxy();
//a*ymax+b=0; a=-getmaxy()/(ymax-ymin)=-1/sc_y; b=-a*ymax=ymax/sc_y;
xx=x/sc_x-x_min/sc_x;
yy=-y/sc_y+y_max/sc_y;
yy1=-ye_0/sc_y+y_max/sc_y;
yy2=-yrk_0/sc_y+y_max/sc_y;
putpixel(xx,yy, 2);
putpixel(xx,yy1,15);
putpixel(xx,yy2,4);
}
getchar();
//closegraph();
}