Решение дифференциальных уравнений по методу Эйлера
Список использованной литературы
1. Д. Мак-Кракен, У. Дорн. Численные методы и программирование на фортране. -М.: Мир,1977.-389,396-408 с.
2. А.А. Самарский. Введение в численные методы. - М.:Наука,1987.-176 с.
3. Алгоритмы вычислительной математики: Лабораторный практикум по курсу «Программирование» для студентов 1 - 2-го курсов всех специальностей БГУИР/А.К. Синицын, А.А. Навроцкий.- Мн.: БГУИР, 2002.- 65-69 с.
4. ГОСТ 2.105-95. Общие требования к текстовым документам.
5. ГОСТ 7.32-91. Система стандартов по информации, библиотечному и издательскому делу. Отчет о НИР. Структура и правила оформления.
Приложение 1. Текст программы.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <conio.h>
#include <dos.h>
#include <graphics.h>
#include <math.h>
#include <bios.h>
//---------------------------------------------------------------------------
void formyl(int p)
{
if(p==1) printf(" 1. C1*y = C2*y + C3*x + C4*x*y");
else if(p==2) printf(" 2. y/(C1-100) = C2*y + C3*x + (C4+x)*y");
else if(p==3) printf(" 3. pow(e,C1)*y = C2*y + C3*cos(x) + (C4+x+y)");
else if(p==4) printf(" 4. C1*sin(x)*y = e*C2*y + C3*arcsin(x) + C4*y/x");
else if(p==5) printf(" 5. C1*y = sin(C2)*y + tg(C3*x) + C4*ln(x)*y");
else if(p==6) printf(" 6. C1*y = y*C2 + C3*sin(x) + C4*cos(x)*y");
else if(p==7) printf(" 7. (C1+C2+C3+C4)*y = C2*y + (C3-x) + lg(C4*x)*y");
else if(p==8) printf(" 8. y/C1 = y/C2 + C3*sin(x) + C4*x*y");
else if(p==9) printf(" 9. sin(C1)*y = C2*y + |C3|*x + x*y/C4");
}
//---------------------------------------------------------------------------
void formyl2(int p,double C1,double C2,double C3,double C4)
{
if(p==1) printf("%.2f*y=%.2f*y+%.2f*x+%.2f*x*y",C1,C2,C3,C4);
else if(p==2) printf("y/(%.2f-100)=%.2f*y+%.2f*x+(%.2f+x)*y",C1,C2,C3,C4);else if(p==3) printf("pow(e,%.2f)*y=%.2f*y+%.2f*cos(x)+(%.2f+x+y)",C1,C2,C3,C4);
else if(p==4) printf("%.2f*sin(x)*y=e*%.2f*y+%.2f*arcsin(x)+%.2f*y/x",C1,C2,C3,C4);
else if(p==5) printf("%.2f*y=sin(%.2f)*y+tg(%.2f*x)+%.2f*ln(x)*y",C1,C2,C3,C4);
else if(p==6) printf("%.2f*y=y*%.2f+%.2f*sin(x)+%.2f*cos(x)*y",C1,C2,C3,C4);
else if(p==7) printf("(%.2f+%f+%.2f+%.2f)*y=%.2f*y+(%.2f-x)+lg(%.2f*x)*y",C1,C2,C3,C4,C2,C3,C4);
else if(p==8) printf("y/%.2f=y/%.2f+%.2f*sin(x)+%.2f*x*y",C1,C2,C3,C4);
else if(p==9) printf("sin(%.2f)*y=%.2f*y+|%.2f|*x+x*y/%.2f",C1,C2,C3,C4);
}
//---------------------------------------------------------------------------
double formyl3(int p,double h,double x,double y,double C1,double C2,double C3,double C4)
{
double Y;
if(p==1) Y=h*(C2*y+C3*x+C4*x*y)/C1+y;
else if(p==2) Y=h*(C1-100)*(y*C2+C3*x+(C4+x)*y)+y;
else if(p==3) Y=h*(C2*y+C3*cos(x)+C4+x+y)/exp(C1)+y;
else if(p==4) Y=h*(exp(1)*C2*y+C3*asin(x)+C4*y/x)/(C1*sin(x))+y;
else if(p==5) Y=h*(sin(C2)*y+tan(C3*x)+C4*log10(x)*y)/C1+y;
else if(p==6) Y=h*(y*C2+C3*sin(x)+C4*cos(x)*y)/C1+y;
else if(p==7) Y=h*(C2*y+(C3-x)+log10(C4*x)*y)/(C1+C2+C3+C4)+y;
else if(p==8) Y=h*(y/C2+C3*sin(x)+C4*x*y)*C1+y;
else if(p==9) Y=h*(C2*y+abs(C3)*x+x*y/C4)/sin(C1)+y;
return(Y);
}
//---------------------------------------------------------------------------
void main(void)
{
int vv=0,vv1=0; // руководит операциями
int N=0,W; // кол промежутков
int i,j,k; // используются во всех "for"
int nom; // номер примера
int st=4,vst=0; // строчка в меню
double C1,C2,C3,C4; // константы
double M; // масштаб
double xtoch,ytoch; // считает y(x) по графику
double h,H; // шаг
double A=0,B=0,ii,jj,kk; // пределы интегрирования
double x[102],y[102]; // главные переменные x,y
//----Подключение графики и проверка-----------------------------------------
int g_driver=9,g_mode=2, g_error;
initgraph(&g_driver,&g_mode,"c:BORLANDgi");
g_error=graphresult();
if(g_error!=grOk)
{
puts("error");
printf(" error=%d, reason=%s ", g_error, grapherrormsg(g_error));
getch();
exit(1);
}
closegraph();
//----Проверка или создание файла--------------------------------------------
FILE *fail;
if((fail = fopen("form.txt", "r")) == NULL)
if((fail = fopen("form.txt", "w")) == NULL)
{
clrscr();
printf("Ошибка при открытии файла");
getch();
exit;
}
fclose(fail);
//----Графическое меню-------------------------------------------------------
menu:
initgraph(&g_driver,&g_mode,"c:BORLANDgi");
cleardevice();
setbkcolor(5);
setfillstyle(7,8); bar(0,0,getmaxx(),getmaxy());
setfillstyle(1,1); bar(200,45,465,145);
setfillstyle(1,15); bar(203,48,462,142);
setfillstyle(1,1); bar(206,51,459,139);
setcolor(11);
settextstyle(7,0,8);
outtextxy(220,40,"Menu");
setcolor(7);
settextstyle(1,0,4);
outtextxy(100,200,"1.Formyla");
outtextxy(100,240,"2.Reshenie");
outtextxy(100,280,"3.Graphic");
outtextxy(100,320," Exit");
settextstyle(2,0,5);
outtextxy(500,440,"beta ver. 101 c");
settextstyle(1,0,4);
izm: // Выбор пункта меню
if(st==1) goto d1;
else if(st==11) goto d11;
else if(st==12) goto d12;
else if(st==2) goto d2;
else if(st==3) goto d3;
else if(st==4) goto d4;
d1:
setcolor(2);
outtextxy(100,200,"1.Formyla");
setcolor(7);
outtextxy(100,240,"2.Reshenie");
outtextxy(100,320," Exit");
if(W==0)
{
W=1;
setfillstyle(7,1);
for(i=0;i<35;i++)
{
bar(280+i*7,190,287+i*7,260);
delay(5);
}
setfillstyle(7,8);
}
outtextxy(300,185,"Enter");
outtextxy(300,215,"Open in fails");
vst=getch();
if(vst==72)
{
st=4;
for(i=0;i<35;i++)
{
bar(530-i*7,190,523-i*7,260);
delay(5);
}
goto izm;
}
else if(vst==80)
{
st=2;
for(i=0;i<35;i++)
{
bar(530-i*7,190,523-i*7,260);
delay(5);
}
goto izm;
}
else if(vst==77) { st=11; goto izm; }
else goto d1;
d11:
setfillstyle(7,1);
bar(280,190,530,260);
setfillstyle(7,8);
setcolor(2);
outtextxy(300,185,"Enter");
setcolor(7);
outtextxy(100,200,"1.Formyla");
outtextxy(100,240,"2.Reshenie");
outtextxy(100,320," Exit");
outtextxy(300,215,"Open in fails");
vst=getch();
if(vst==72) { st=12; goto izm; }
else if(vst==80) { st=12; goto izm; }
else if(vst==75) { st=1; goto izm; }
else if(vst==13) { vv=1; vv1=1; closegraph(); goto resh; }
else goto d11;
d12:
setfillstyle(7,1);
bar(280,190,530,260);
setfillstyle(7,8);
setcolor(2);
outtextxy(300,215,"Open in fails");
setcolor(7);
outtextxy(100,200,"1.Formyla");
outtextxy(100,240,"2.Reshenie");
outtextxy(100,320," Exit");
outtextxy(300,185,"Enter");
vst=getch();
if(vst==72) { st=11; goto izm; }
else if(vst==80) { st=11; goto izm; }
else if(vst==75) { st=1; goto izm; }
else if(vst==13) { vv=1; vv1=2; closegraph(); goto resh; }
else goto d12;
d2:
setfillstyle(1,5);
setfillstyle(7,8);
bar(280,160,600,400);
setcolor(2);
outtextxy(100,240,"2.Reshenie");
setcolor(7);
outtextxy(100,200,"1.Formyla");
outtextxy(100,280,"3.Graphic");
vst=getch();
if(vst==13) { vv=2; closegraph(); goto resh; }
else if(vst==72) { st=1; W=0; goto izm; }
else if(vst==80) { st=3; goto izm; }
else goto d2;
d3:
setcolor(2);
outtextxy(100,280,"3.Graphic");
setcolor(7);
outtextxy(100,240,"2.Reshenie");
outtextxy(100,320," Exit");
vst=getch();
if(vst==13) { vv=3; closegraph(); goto resh; }
else if(vst==72) { st=2; goto izm; }
else if(vst==80) { st=4; goto izm; }
else goto d3;
d4:
setfillstyle(1,5);
setfillstyle(7,8);
bar(280,160,600,400);
setcolor(2);
outtextxy(100,320," Exit");
setcolor(7);
outtextxy(100,200,"1.Formyla");
outtextxy(100,240,"2.Reshenie");
outtextxy(100,280,"3.Graphic");
vst=getch();
if(vst==13) { vv=4; closegraph(); goto resh; }
else if(vst==72) { st=3; goto izm; }
else if(vst==80) { st=1; W=0; goto izm; }
else goto d4;
//----Ввод данных------------------------------------------------------------
resh:
if(vv==1)
{
if(vv1==1)
{
C1=0;C2=0;C3=0;C4=0;
clrscr();
printf(" Выбор формулы решения ");
for(i=1;i<10;i++)
formyl(i);
printf(" нажмите любой номер: ");
while(2)
{
nom=getch()-48;
if((nom>0)&(nom<10))
break;
}
while(3)
{
clrscr();
printf(" по окончании ввода констант нажать - Esc");
printf(" выбронная формула имеет вид ");
formyl(nom);
printf(" ввод констант в формулу дифференциального уравнения: ");
formyl2(nom,C1,C2,C3,C4);
printf(" C");
i=getch();
if(i==49)
{
printf("1 = ");
scanf("%lf",&C1);
}
else if(i==50)
{
printf("2 = ");
scanf("%lf",&C2);
}
else if(i==51)
{
printf("3 = ");
scanf("%lf",&C3);
}
else if(i==52)
{
printf("4 = ");
scanf("%lf",&C4);
}
else if(i==27)
{
break;
}
}
AB:
clrscr();
printf(" выбронная формула имеет вид ");
formyl(nom);
printf(" ввод констант в формулу дифференциального уравнения: ");
formyl2(nom,C1,C2,C3,C4);
printf(" уточните границы изменения Х, от A= ");
scanf("%lf",&A);
clrscr();
printf(" выбронная формула имеет вид ");
formyl(nom);
printf(" ввод констант в формулу дифференциального уравнения: ");
formyl2(nom,C1,C2,C3,C4);
printf(" уточните границы изменения Х, от A= %.2f до B= ",A);
scanf("%lf",&B);
if(B<=A)
{
printf(" неправильно заданны границы!!! ");
delay(1000);
goto AB;
}
clrscr();
printf(" выбронная формула имеет вид ");
formyl(nom);
printf(" ввод констант в формулу дифференциального уравнения: ");
formyl2(nom,C1,C2,C3,C4);
printf(" уточните границы изменения Х, от A= %.2f до B= %.2f",A,B);
printf(" задайте начальное значение Y(%.2f) = ",A);
scanf("%lf",&y[0]);
printf(" сохранить данные в файле (Да/Нет) ");
while(11)
{
i=getch();
if((i==76)||(i==108)||(i==132)||(i==164))
{
printf("Да");
delay(300);
fopen("form.txt", "w");
fwrite(&nom,1,2,fail);
fwrite(&C1,1,8,fail);
fwrite(&C2,1,8,fail);
fwrite(&C3,1,8,fail);
fwrite(&C4,1,8,fail);
fwrite(&A,1,8,fail);
fwrite(&B,1,8,fail);
fwrite(&y[0],1,8,fail);
fclose(fail);
break;
}
else if((i==89)||(i==121)||(i==141)||(i==173))
{
printf("Нет");
delay(300);
break;
}
}
goto menu;
}
if(vv1==2)
{
fopen("form.txt", "r");
fread(&nom,1,2,fail);
fread(&C1,1,8,fail);
fread(&C2,1,8,fail);
fread(&C3,1,8,fail);
fread(&C4,1,8,fail);
fread(&A,1,8,fail);
fread(&B,1,8,fail);
fread(&y[0],1,8,fail);
fclose(fail);
clrscr();
printf(" выбранная формула имеет вид ");
formyl(nom);
printf(" ввод констант в формулу дифференциального уравнения: ");
formyl2(nom,C1,C2,C3,C4);
printf(" уточните границы изменения Х, от A= %.2f до B= %.2f",A,B);
printf(" задайте начальное значение Y(%.2f) = %.2f",A,y[0]);
printf(" данные были прочитаны из ранее сохраненного файла");
getch();
goto menu;
}
}
//----Решение примера--------------------------------------------------------
else if(vv==2)
{
if((A==0)&(B==0))
{
clrscr();
printf(" Сперва необходимо записать пример!");
delay(2000);
goto menu;
}
clrscr();
printf(" для решения этой задачи нужно: ");
printf(" задать колличество отрезков, на которое расделится промежуток AB=[%.2f,%.2f] N= ",A,B);
scanf("%d",&N);
ii=B; jj=A;
h=(B-A)/N;
N++;
//-------------------------------------
for(i=0;i<N;i++)
{
x[i]=A+h*i;
y[i+1]=formyl3(nom,h,x[i],y[i],C1,C2,C3,C4);
}
clrscr();
printf(" Решение задачи: ");
printf("
г========T===================T========================¬ ");
printf(" ¦ i ¦ x=h*i ¦ y(i) ¦ ");
printf(" ¦========+===================+========================¦ ");
for(i=0;i<N;i++)
printf(" ¦ %3d ¦ %11lf ¦ %14f ¦ ",i+1,x[i],y[i]);
printf(" L========¦===================¦========================- ");
getch();
}
//----График-----------------------------------------------------------------
else if(vv==3)
{
if(N==0)
{
clrscr();
printf(" Сперва необходимо решить пример!");
delay(2000);
goto menu;
}
initgraph(&g_driver,&g_mode,"c:BORLANDgi");
cleardevice();
setbkcolor(0);
setcolor(11);
settextstyle(4,0,5);
outtextxy(100,280,"Open graphic...");
setfillstyle(1,7);
bar(218,238,422,252);
setfillstyle(1,5);
bar(220,240,420,250);
setcolor(9);
for(i=0;i<200;i=i+2)
{
line(221+i,241,221+i,249);
delay(10);
}
M=5;
while(31)
{
// Фон
clrscr();
cleardevice();
setbkcolor(0);
setcolor(10);
settextstyle(5,0,4);
// Вывод данных по графику
printf(" выход в меню - Esc");
printf(" график дифференциального уравнения - ");
formyl2(nom,C1,C2,C3,C4);
printf(" построенный с шагом h = %f",h);
printf(" на промежутке от %.2f до %.2f",A,B);
printf(" выполненный в масштабе %.4f: 1",M);
printf(" найти решение в точке - Enter");
printf(" изменение масштаба увеличить - +");
printf(" уменьшить - -");
// обводка слов
setcolor(3);
line(455,0,630,0);
line(455,15,630,15);
line(455,0,455,15);
line(630,0,630,15);
setcolor(9);
// Оси
setcolor(10);
line(50,250,600,250);
line(580,245,600,250);
line(580,255,600,250);
line(325,50,325,440);
line(325,50,330,70);
line(325,50,320,70);
setcolor(15);
settextstyle(1,0,1);
outtextxy(600,260,"X");
outtextxy(345,50,"Y");
outtextxy(300,260,"O");
// единичные насечки
if(M>4)
for(i=1;i<70;i++)
{
if(i*M<=255)
{
line(325+M*i,248,325+M*i,252);
line(325-M*i,248,325-M*i,252);
}
if(i*M<=170)
{
line(323,250-M*i,327,250-M*i);
line(323,250+M*i,327,250+M*i);
}
}
// Прорисовка графика
setcolor(9);
for(i=1;i<N;i++)
{
if((M*y[i]>-1000)&(M*y[i]<1000))
{
line(325+M*((i-1)*h+A),250-M*y[i-1],325+M*(i*h+A),250-M*y[i]);
delay(10);
}
}
k=getch();
if(k==43) // масштаб +
{
if (M>=1) M++;
else if((M<1) &&(M>=0.1)) M=M+0.1;
else if((M<0.1) &&(M>=0.01)) M=M+0.01;
else if((M<0.01) &&(M>=0.001)) M=M+0.001;
else if((M<0.001) &&(M>=0.0001)) M=M+0.0001;
else if((M<0.0001) &&(M>=0)) M=M+0.00001;
}
else if(k==45) // масштаб -
{
if (M>1.000001) M--;
else if((M<=1) &&(M>0.100001)) M=M-0.1;
else if((M<=0.100001)&&(M>0.010001)) M=M-0.01;
else if((M<=0.010001)&&(M>0.001001)) M=M-0.001;
else if((M<=0.001001)&&(M>0.000101)) M=M-0.0001;
else if((M<=0.000100)&&(M>0.000010)) M=M-0.00001;
}
else if(k==13) // Нахождение y по значению x
{
printf(" введите значение абсциссы x= ");
scanf("%lf",&xtoch);
if((A<=xtoch)&&(xtoch<=B))
{
setcolor(11);
for(i=1;i<N;i++)
{
if((x[i-1]<=xtoch)&&(xtoch<=x[i]))
{
H=xtoch-x[i-1];
ytoch=formyl3(nom,H,x[i-1],y[i-1],C1,C2,C3,C4);
printf("ордината точки равна y= %lf",ytoch);
line(325+M*xtoch,250,325+M*xtoch,250-M*ytoch);
break;
}
}
setcolor(9);
}
else
printf("введена абсцисса, не пренадлежащая промежутку решения");
getch();
}
else if(k==27)
break;
}
closegraph();
}
//----Выход из программы-----------------------------------------------------
else if(vv==4)
{
exit(2);
}
//----Ввод программы в бесконечный цикл--------------------------------------
goto menu;
}
Приложение 2. Результаты работы программы
Приведем расчет дифференциального уравнения первого, второго и третьего порядка методом Эйлера
1. Пусть дано дифференциальное уравнение первого порядка: y/=2x-y
Требуется найти решение на отрезке [0,1] c шагом h=(1-0)/5=0,2
Начальные условия: у0=1;
Пользуясь рекурентными формулами (4), находим:
1). x1=0,2; х1/2=0,1; y(x1)=y(x0)+б0h; y(x1/2)=y(x0)+f(x0,y0)h/2;
f(x0,y0)=20-1=-1
y(x1/2)=1-10,1=0,9
б0=20,1-0,9=-0,7
y1=1-0,10,2=0,86
2). y(x2)=y(x1)+б1h; x2=0,2+0,2=0,4; x1+1/2=x1+h/2=0,2+0,1=0,3
y(x1+1/2)=y(x1)+f(x1,y(x1))h/2
f(x1,y1)=20,2-0,86=-0,46
y(x1+1/2)=0,86-0,460,1=0,814
б1=2*0,3-0,814=-0,214
y2=0,86-0,214*0,2=0,8172
3). x3=0,4+0,2=0,6; x2+1/2=x2+h/2=0,4+0,1=0,5
f(x2,y2)=2*0,4-0,8172=-0,0172
y2+1/2=0,8172-0,0172*0,1=0,81548
б2=2*0,5-0,81548=0,18452
y3=0,8172+0,18452*0,2=0,854104
4).x4=0,8; x3+1/2=x3+h/2=0,6+0,1=0,7
f(x3,y3)=2*0,6-0,854104=0,345896
y3+1/2=0,854104+0,345896*0,1=0,8886936
б3=2*0,7-0,89=0,5113064
y4=0,854104+0,5113064*0,2=0,95636528
5).x5=1; x4+1/2=0,8+0,1=0,9
f(x4,y4)=2*0,8-0,956=0,64363472
y4+1/2=0,956+0,643*0,1=1,020728752;
б4=2*0,9-1,02=0,779271248
y5=0,956+0,7792*0,2=1,11221953