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
правее узловых точек.
При
epsilon = 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 Pascal, то должен получиться примерно такой
текст:
Uses crt;
function dy(x,y:real):real;
begin dy:= y+exp(x); end; var
x0,y0,x1,y,h,x,dydx :real;
i,n
:integer;
begin
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. |