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

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

"Самиздат": [Регистрация] [Найти] [Рейтинги] [Обсуждения] [Новинки] [Обзоры] [Помощь|Техвопросы]
Ссылки:


 Ваша оценка:
  • Аннотация:
    Неплохой новый метод для численного решения обычных дифференциальных уравнений


   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 при котором интеграл будет вычислен более точно.
   Здесь я не буду описывать всю процедуру поиска решения, который я когда-то прошёл (хотя, это было бы многим интересно), а расскажу самую суть этого метода.
   Итак, нас с вами надо вычислить интеграл:
  
   0x01 graphic
   рис.1
  
   Будем вычислять интеграл методом левых прямоугольников:
  
   n = 4
   h = (b-a)/n
   int1 = 0
   for i=0 to n-1
   int1 = int1 +h*fnF(x0+i*h)
   next i
  
   То есть, функцию мы вычисляем в левых точках:
   0x08 graphic
  
  
   Рис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 эта функция графически будет выглядеть как прямая линия:
  
   0x01 graphic
   рис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
  
   Кто-то сейчас наверно задаёт вопрос: "И что это даёт?". Вопрос хороший.
   Теперь вычислим тот же интеграл другим способом:
  
   0x01 graphic
  
   рис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, а интеграл уравнения (**), который выглядит так :
  
   0x01 graphic
   рис5
  
  
   Левая часть интеграла (***) как раз и равна fi1 или fi0 , взависимости от того, чему равен эпсилон. Запись y(x1) есть тоже самое, что и запись y1 .
  
   Несколько сложнее вычислять betta1 и betta0. Для этого надо (***) переписать так:
  
   0x01 graphic
  
   рис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 года

 Ваша оценка:

Связаться с программистом сайта.

Новые книги авторов СИ, вышедшие из печати:
Э.Бланк "Пленница чужого мира" О.Копылова "Невеста звездного принца" А.Позин "Меч Тамерлана.Крестьянский сын,дворянская дочь"

Как попасть в этoт список
Сайт - "Художники" .. || .. Доска об'явлений "Книги"