Естественно,
чтобы вычислить интегралы fi1 и
fi0 , надо эту программу немного переделать
:
uses crt ;
function dy(x,y:real):real; begin dy:= y+exp(x); end;
function dy2(x,y:real):real;begin dy2:=dy(x,y)+exp(x); end;
var
i,n :integer;
x0,x1,y0,y1, h,x,dydx,
y,epsi,f1, fi0,fi1,be0,be1,z :real;
begin clrscr;
x0:= 1 ;
x1:= 2.4 ;
y0:= 2.71828;
writeln('x0=',x0:5:2,' y0= ',y0:7:5);
n:= 4;
writeln('N=',n);
h:=(x1-x0)/n;
epsi:= 0.0;
y:= y0+epsi*h*dy(x0,y0);
for i:=1 to n do
begin
x:=x0+(i-1)*h+epsi*h;
dydx:= dy(x,y);
y:=y+h*dydx;
end;
writeln('epsi=',epsi:5:2,' y1= ',y:12:6);
fi0:= y-y0;
writeln(' fi0= ',fi0:12:6);
writeln;
epsi:= 0.999;
y:= y0+epsi*h*dy(x0,y0);
for i:=1 to n do
begin
x:=x0+(i-1)*h+epsi*h;
dydx:= dy(x,y);
y:=y+h*dydx;
end;
writeln('epsi=',epsi:5:2,' y1= ',y:12:6);
fi1:= y-y0;
writeln(' fi1= ',fi1:12:6);
if readkey=#3 then;end.
Связано это с
тем, что нас интересует не точка y1, а интеграл
уравнения (**), который выглядит так :
рис5
Левая часть
интеграла (***) как раз и равна fi1 или fi0 , взависимости от того, чему равен эпсилон. Запись y(x1)
есть тоже
самое, что и запись y1
.
Несколько сложнее вычислять betta1 и betta0. Для этого надо (***) переписать так:
рис6
После этого уравнение (****) надо решать либо численно
относительно y1 , либо аналитически, если это
возможно.
А затем, после этого, уже можно
вычислить:
betta0 = y1 – y0
или
betta1 = y1 – y0
взависимости
от того, чему равен эпсилон.
Вся программа
решения уравнения (**) выглядит так :
uses crt ;
function f(x:real):real; begin
f:=x*exp(x); end;
function dy(x,y:real):real; begin dy:= y+exp(x); end;
function dy2(x,y:real):real;begin dy2:=dy(x,y)+exp(x); end;
var
i,n :integer;
x0,x1,y0,y1, h,x,dydx,
y,epsi,f1, fi0,fi1,be0,be1,z :real;
begin clrscr;
textcolor(3);
x0:= 1 ;
x1:= 2.4 ;
y0:= f(x0); writeln('x0=',x0:5:2,' y0= ',y0:7:5);
y1:= f(x1); writeln('Точно: y1= ',y1:12:5);
n:= 4;
writeln('N=',n);
h:=(x1-x0)/n;
y:= y0;
for i:=1 to n do
begin
x:=x0+(i-1)*h;
dydx:= dy(x,y);
y:=y+h*dydx;
end;
textcolor(14);
writeln(' Метод Эйлера
даёт: y1= ',y:12:6);
textcolor(3);
writeln('=============================');
epsi:= 0.0;
y:= y0+epsi*h*dy(x0,y0);
for i:=1 to n do
begin
x:=x0+(i-1)*h+epsi*h;
dydx:= dy(x,y);
y:=y+h*dydx;
end;
writeln('epsi=',epsi:5:2,' y1= ',y:12:6);
fi0:=y-y0;
writeln(' fi0= ',fi0:12:6);
writeln;
epsi:= 0.999;
y:= y0+epsi*h*dy(x0,y0);
for i:=1 to n do
begin
x:=x0+(i-1)*h+epsi*h;
dydx:= dy(x,y);
y:=y+h*dydx;
end;
writeln('epsi=',epsi:5:2,' y1= ',y:12:6);
fi1:=y-y0;
writeln(' fi1= ',fi1:12:6);
writeln('=============================');
y:=y0;
z:=0;
for i:=1 to n do
begin
x:=x0+(i-1)*h;
z:=z +h*x*dy2(x,y);
y:=y+h*dy(x,y);
end;
y1:=
(y0-x0*dy(x0,y0)-z+x1*exp(x1) )/(1-x1);
writeln('epsi=0 y1=',y1:9:4);
be0:=y1-y0;
writeln(' be0=',be0:9:4);
writeln;
epsi:=0.999;
y:=y0+epsi*h*dy(x0,y0);
z:=0;
for i:=1 to n do
begin
x:=x0+(i-1)*h +epsi*h;
z:=z +h*x*dy2(x,y);
y:=y+h*dy(x,y);
end;
y1:=
(y0-x0*dy(x0,y0)-z+x1*exp(x1) )/(1-x1);
writeln('epsi=1 y1=',y1:9:4);
be1:= y1-y0;
writeln(' be1=',be1:9:4);
writeln('=============================');
z:= fi0*(be1-be0)-be0*(fi1-fi0);
z:=z/(be1-be0 -(fi1-fi0) );
textcolor(14);
writeln(' Эпсилон-метод: y1 =',z+y0:9:4);
if
readkey=#3 then;end.
Возможно,
кому-то этот метод не показался слишком красивым. А на самом деле – на практике
– все методы хороши. Иной раз, там, где не очень удобен самый хороший, самый
"дорогой" способ, другой метод, неказистый и слишком простой, отлично
решает поставленную перед ним задачу…
Но, честно
говоря, предложенный мной метод не является более точным, чем метод
Рунге-Кутта.
Павел Сапунов
1 декабря 2007 года |