1. Эпсилон-метод
Есть такой метод численного интегрирования: сначала готовится обычный способ численного интегрирования методом прямоугольников:
'' листинг 1
n=8
h=(x1-x0)/n
Integral = 0
For i=0 to n-1
Integral = integral +h*fnF(x0+i*h)
Next i
'' язык программирования - Бейсик
'' h - шаг численного интегрирования
а затем подынтегральная функция вычисляется не в узлах, а в точках между ними:
'' листинг 2
n=8
h=(x1-x0)/n
Integral = 0
For i=0 to n-1
Integral = integral +h*fnF(x0+i*h+0.5*h)
Next i
Этот метод даёт неплохой результат: например, при интегрировании функции fnF(x)=x^2 от 1 до 3 ( x0=1 ; x1=3 ) при n=5 мы получаем такой же результат, как и обычным методом прямоугольников при n=300 (неплохо для такого простого метода!). И поэтому я, осознав однажды этот факт, решил усовершенствовать этот итак не плохой метод. И у меня всё получилось.
1.1 Численное интегрирование
Как же мне удалось в корне изменить данный метод? Во-первых, вместо строчки
Integral = integral +h*fnF(x0+i*h+0.5*h)
я написал такую:
Integral = integral +h*fnF(x0+i*h+epsilon*h)
(к сожалению греческая буква эпсилон не на каждом сайте печатается)
Заменить число на букву - это в математике большое дело... Не верите? Хорошо, читайте дальше - я вам докажу.
Наверно, Вы уже догадались, что я решил найти вместо 0.5 такое число epsilon при котором интеграл будет вычислен более точно.
Здесь я не буду описывать всю процедуру поиска решения, который я когда-то прошёл (хотя, это было бы многим интересно), а расскажу самую суть этого метода.
Итак, нас с вами надо вычислить интеграл:
рис.1
Будем вычислять интеграл методом левых прямоугольников:
n = 4
h = (b-a)/n
int1 = 0
for i=0 to n-1
int1 = int1 +h*fnF(x0+i*h)
next i
То есть, функцию мы вычисляем в левых точках:
Рис2
Пусть, для начала, эпсилон будет постоянной величиной:
epsilon=const ( 0<= epsilon <=1 )
Если эпсилон равно нулю, то это будет соответствовать листингу 1, если эпсилон=0.5, то это будет соответствовать листингу 2. А если эпсилон будет равен 0.3, то интеграл будет вычисляться в точках, смещенных на 0.3h правее узловых точек.
При эпсилон=const алгоритм интегрирования будет таким:
'' листинг 3
n=8
h=(x1-x0)/n
Integral = 0
For i=0 to n-1
Integral = integral +h*fnF(x0+i*h+ ?*h)
Next i
А если эпсилон изменять плавно от нуля до единицы? Что тогда будет? Тогда у нас будет иметь место функция Фи от эпсилон:
fi = fi( epsilon )
При большом n эта функция графически будет выглядеть как прямая линия:
рис3
Увидеть эту линию можно, набрав небольшую программку:
'' listing 4
screen 9
def fnF(x)=x^3
x00=100
y00=200
x0=1
x1=3
n=10
h=(x1-x0)/n
xPoint=0
for epsilon=0 to 1 step 0.01
integral=0
for i=0 to n-1
integral=integral+fnF(x0+i*h+epsilon*h)*h
next i
pset(x00+xPoint, y00-integral*4),2
xPoint=xPoint+1
next epsilon
end
Кто-то сейчас наверно задаёт вопрос: "И что это даёт?". Вопрос хороший.
Теперь вычислим тот же интеграл другим способом:
рис4
При аналитическом интегрировании формула (2) даёт такой же результат, что и формула (1). Но при численном интегрировании формула (2) - это совсем иная вещь!
При численном интегрировании формулы (1) и (2) обладают замечательным свойством: графики этих двух функций - fi(epsilon ) и betta(epsilon) - являются почти прямыми линиями (при большом n ), причём они всегда пересекаются!
Точку их пересечения назовём epsilon* . И сразу скажем, что эту точку можно вычислить. Причём, очень просто:
epsilon* = (betta0 - fi0)/( fi1 - fi0 - betta1 + betta0 ) (3)
где:
betta0 - интеграл по формуле (2) при epsilon=0
betta1 - интеграл по формуле (2) при epsilon =1
fi0 - интеграл по формуле (1) при epsilon =0
fi1 - интеграл по формуле (1) при epsilon =1
Но если точка epsilon* - это тот самый замечательный эпсилон, в котором интеграл будет наиболее точным, значит можно вычислить интеграл как минимум двумя способами: во-первых, можно использовать программу листинга 3, во-вторых, можно просто сделать так:
integral* = ( G fi0 - betta0 )/( G - 1) (4)
где
G = ( -b F'(b) + a F'(a) )/( F(b) - F(a) )
Где F'(b) - производная функции F(х) в точке b
F'(a) - производная функции F(х) в точке a
Самые догадливые из читателей наверно уже догадались, где именно в этом методе есть слабое место. Да, вы правы, при малых значениях n (n = 1,2,3...) функция fi(epsilon ) ( также как и betta(epsilon ) ) не является похожей на прямую линию. Кстати, у любого из вас есть шанс усовершенствовать ЭПСИЛОН-МЕТОД: ведь подобрать функцию betta(epsilon ) при малых n не так уж сложно... Дерзайте! Удачи!
Но когда будете получать премию за самый лучший численный метод, на забудьте сказать пару добрых слов и про меня, скромного любителя численных методов...
1.2 Численное решение дифференциальных уравнений
Для того, чтобы решить дифференциальное уравнение с помощью эпсилон-метода, надо потрудиться. Потому что с помощью этого метода решать "дифуры" не очень легко.
Скажем надо решить диф.ур-ие:
y'x = y+exp(x)
в общем виде
y'x = F(x,y) (**)
начальные условия, например, такие:
x0 = 1
y0 = 2.71828
Перед тем, как решать, надо договориться, какую формулу мы будем использовать: (3) или (4).
Используя своё особое положение (я - автор), я предлагаю использовать (4). Итак, используя ос.пол., начинаю использовать всё остальное...
А именно: имею право использовать новую формулу:
integral* = ( fi0 (betta1 - betta0 ) - betta0 ( fi1 - betta0) ) / ( (betta1 - betta0 ) - (fi1 - fi0) ) (5)
Она самая удобная.
Тут я, прийдя к выводу, что язык Turbo Basic очень любит делать ошибки в вычислениях, я решил перейти к языку Turbo Pascal. Если решать уравнение (**) методом Эйлера с помощью Turbo Pascal, то должен получиться примерно такой текст:
Uses crt;
function dy(x,y:real):real;
begin dy:= y+exp(x); end;
x0:=1; y0:=2.71828;
x1:=2.4;
n:=4;
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;
writeln( 'Методом Эйлера: y1=', y:9:4);
readln;end.
Естественно, чтобы вычислить интегралы 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 года