Сапунов Павел. Численные методы. Метод "Эпсилон" для решения дифференциальных уравнений

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* = (betta₀- fi₀)/( fi₁ - fi₀ - betta₁ + betta₀ ) (3)

где:

betta₀- интеграл по формуле (2) при epsilon=0

betta₁- интеграл по формуле (2) при epsilon =1

fi₀- интеграл по формуле (1) при epsilon =0

fi₁- интеграл по формуле (1) при epsilon =1

Но если точка epsilon* - это тот самый замечательный эпсилон, в котором интеграл будет наиболее точным, значит можно вычислить интеграл как минимум двумя способами: во-первых, можно использовать программу листинга 3, во-вторых, можно просто сделать так:

integral* = ( G fi₀- betta₀)/( 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) (**)

начальные условия, например, такие:

x₀ = 1

y₀ = 2.71828

Перед тем, как решать, надо договориться, какую формулу мы будем использовать: (3) или (4).

Используя своё особое положение (я - автор), я предлагаю использовать (4). Итак, используя ос.пол., начинаю использовать всё остальное...

А именно: имею право использовать новую формулу:

integral* = ( fi₀(betta₁- betta₀) - betta₀( fi1 - betta₀) ) / ( (betta₁- betta₀) - (fi₁ - fi₀) ) (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.

Естественно, чтобы вычислить интегралы fi₁и fi₀ , надо эту программу немного переделать :

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) есть тоже самое, что и запись y₁ .

Несколько сложнее вычислять betta_{1 и}betta_0.Для этого надо (***) переписать так:

рис6

После этого уравнение (****) надо решать либо численно относительно y₁ , либо аналитически, если это возможно.

А затем, после этого, уже можно вычислить:

betta₀= y₁ - y₀

или

betta₁= y₁ - y₀

взависимости от того, чему равен эпсилон.

Вся программа решения уравнения (**) выглядит так :

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 года

Комментарии: 2, последний от 31/12/2011. © Copyright Сапунов Павел (pav.sapunov@ya.ru) Размещен: 24/12/2011, изменен: 24/12/2011. 26k. Статистика. Статья: Естествознание Математические миниатюры Иллюстрации/приложения: 6 шт. Скачать FB2		Ваша оценка:
Аннотация: Неплохой новый метод для численного решения обычных дифференциальных уравнений

Сапунов Павел : другие произведения.

Численные методы. Метод "Эпсилон" для решения дифференциальных уравнений