Численное интегрирование , [Pascal] Решаем интеграл разными способами

[e-moe]

subj in attach (CP-866)

Че скажите?

09.05.05: Исправил ошибку в методе Симпсона...

Сообщение отредактировано: Jin X - 17.11.16, 15:19

Прикреплённый файлLAB3_.zip (1.51 Кбайт, скачиваний: 318)

[Romtek]

Теория - это, конечно, хорошо. Но всё же я рекомендую составлять поэтапное решение поставленных задач, опираясь на источники теории, дабы научить народ переводить математические формулы в программный код.
Важна "история" развития программы: с чего начинать (типы, процедуры, циклы), методы и пути решения, и т.д.

[e-moe]

ok, буду че-то думать

[e-moe]

Долго думал куда это запостить: то ли в "Маразм прогамееров", то ли в топик, где обсуждались самые "корявые исходники"
Решил запостить сюда, ведь это все работает(пусть не так оптимально как хотелось бы, но зато мало строк )....
Тем более

Цитата Romtek @ 02.05.05, 22:08

методы и пути решения, и т.д.

имхо, если грамотно раписать как это "чудо" работает, то это будет полезнее чем просто расписать как программить формулы
Подробнее см. комменты в конце кода...

З.Ы. че скажите?

{
  Writeln('Labinskiy Nikolay aka e-moe (c) 2005');
  Writeln('See you later on forum.sources.ru ;)');
}
{$N+,E+,G+}
const
  a = -1.1;
  b = 0.9;
 
type
  Tfunc = function(x: extended):extended;
 
function func(x: extended): extended; far;
begin
  func:=0.9*cos(-1.77*x)+2.65*exp(-1.57*x)-11.48*x*x*x*x+1.29*x-10.49;
end;
 
function Int1(f: TFunc; a,b: extended; n: byte): extended;
{ Правило трапеций }
var
  res: double;   { времееная переменная, будет хранить результат }
  {если extended, то надо обнулить отдельно, а не в цикле :( TP7.1}
  i: byte;       { счетцик цикла }
begin
  { на 0 делить нельзя! }
  inc(n,byte(n=0));
  for i:=0 to n do
    res:=res*byte(i<>0)+(byte((i>0)and(i<n))+1)*f(a+i*(b-a)/n)*(b-a)/n*0.5;
      { интеграл - площадь фигуры, а следовательно не может
        быть отрицательной, поэтому используем абсолютное значение
        функции (или по-русски модуль) }
  Int1:=res;
end;
 
function Int2(f: TFunc; a,b: extended; n,m: byte): extended;
begin
  { n<>m && n*m <> 0 }
  inc(n,byte(n=0));
  inc(m,byte(m=0)*(byte(n=1)+1)-byte((m=255) and (n=255))*2+byte(m=n));
  { Для подготовки к вычислению надо в 2 раза больше строк чем
    на само вычисление... Странно, не так ли? :) }
  Int2:=Int1(f,a,b,n)+(Int1(f,a,b,n)-Int1(f,a,b,m))/((n*n)/(m*m)-1)
end;
 
function Int3(f: TFunc; a,b: extended; n: byte): extended;
var
  res: double; {если extended, то надо обнулить отдельно, а не в цикле :( TP7.1}
  i: byte;
begin
  inc(n,byte(n and 1)+byte(n=0)*2-byte(n=255)*2);  { (n mod 2=0) && n<>0 }
  for i:=0 to n do
    res:=res*byte(n<>0)+(1 shl (byte((i>0)and(i<n))+byte(odd(i)and(i>0)and(i<n))))*f(a+i*(b-a)/n)*(b-a)/n/3;
  Int3:=res;
end;
 
function Int4(f: TFunc; a,b: extended; n: byte): extended;
const
  AiMi : array[2..6,1..12] of extended =
         ((0.5773502692,1,-0.5773502692,1,0,0,0,0,0,0,0,0),
         (0.7745966692,0.5555555556,-0.7745966692,0.5555555556,
           0,0.8888888889,0,0,0,0,0,0),
         (0.8611363116,0.3478548451,-0.8611363116,0.3478548451,
           0.3399810436,0.6521451549,-0.3399810436,0.6521451549,0,0,0,0),
         (0.9061798459,0.2369268851,-0.9061798459,0.2369268851,0.5384693101,
           0.4786286705,-0.5384693101,0.4786286705,0.0,0.5688888889,0,0),
         (0.9324695142,0.1713244924,-0.9324695142,0.1713244924,
           0.6612093865,0.3607615730,-0.6612093865,0.3607615730,
           0.2386191861,0.4679139346,-0.2386191861,0.4679139346));
var
  i:byte;
  res: double; {если extended, то надо обнулить отдельно, а не в цикле :( TP7.1}
begin
  { n>=2 && n <=6 }
  inc(n,byte(n<2)*(2-n)-byte(n>6)*(n-6));
  for i:=1 to n do
    res:=res*byte(i<>1)+AiMi[n,i*2]*((b-a)*f((b-a)*AiMi[n,i*2-1]/2+(b+a)/2)/2);
  Int4:=res;
end;
 
begin
  Writeln('     Численное интегрирование:',#10);
  Writeln('f(x) = 0.9*cos(-1.77x) + 2.65*e^(-1.57x) - 11.48*x^4 + 1.29x - 10.49');
  Writeln('Int(-1.1, 0.9)[f(x)]',#10);
 
  Writeln('Правило трапеций:');
  writeln('n=2      I=',Int1(func,a,b,2):5:5);
  writeln('n=8      I=',Int1(func,a,b,8):5:5);
  writeln('n=32     I=',Int1(func,a,b,32):5:5);
  writeln('n=64     I=',Int1(func,a,b,64):5:5);
 
  Writeln('Экстраполяционный переход к пределу::');
  writeln('n=2,8    I=',Int2(func,a,b,2,8):5:5);
  writeln('n=8,16   I=',Int2(func,a,b,8,16):5:5);
  writeln('n=16,32  I=',Int2(func,a,b,16,32):5:5);
 
  Writeln('Правило Симпсона:');
  writeln('n=2      I=',Int3(func,a,b,2):5:5);
  writeln('n=8      I=',Int3(func,a,b,8):5:5);
  writeln('n=32     I=',Int3(func,a,b,32):5:5);
 
  Writeln('Метод Гаусса:');
  writeln('n=2      I=',Int4(func,a,b,2):5:5);
  writeln('n=3      I=',Int4(func,a,b,3):5:5);
  writeln('n=4      I=',Int4(func,a,b,4):5:5);
  writeln('n=5      I=',Int4(func,a,b,5):5:5);
  writeln('n=6      I=',Int4(func,a,b,6):5:5);
 
end.
З.Ы.
     не знаю что на меня нашло :), но я решил макисмально
     сократить число строк во всех функциях, вот что из
     всего этого вышло... :)
 
З.З.Ы.
      параметры функций:
        f - ссылка на подинтегральную функцию
        a,b - координаты начальной и конечной точнек интегрирования
        n[,m] - чило разбиений.
 
     в функции можно передавать любые n и m: если нужно,
     то прога сама их подправит так, как нужно для
     конкретного метода;)
 
З.З.З.Ы.
     как это работает? а кто его знает... я уже и
     сам не пойму :) а если серьезно, то могу помочь
     разобраться что тут к чему ;)

Че-то мне кажется, что более короткой реализации этих процедур еще не делали

Сообщение отредактировано: e-moe - 15.06.06, 17:58

[Romtek]

e-moe, а теперь подробнее про каждый метод:
Что, как и почему.
Не знаю на кого рассчитан сей код, но я его с трудом понимаю.

[e-moe]

Цитата Romtek @ 07.05.05, 20:06

e-moe, а теперь подробнее про каждый метод:

Угу, я это и хочу расписать )
Только вот не знаю, может сделать отдельную тему где писать все такие "штучки", а тут все-таки написать нормальный исходник по интегрированию?

Цитата Romtek @ 07.05.05, 20:06

...но я его с трудом понимаю.

Эт точно, если не расписать как это работает, то я и сам скоро забуду....

[Romtek]

Цитата e-moe @ 09.05.05, 09:39

может сделать отдельную тему где писать все такие "штучки", а тут все-таки написать нормальный исходник по интегрированию?

Зачем плодить лишние темы? "Штучки" можно вполне и здесь обсудить (это и есть цель темы).
Как только появится что-то оформленное и обсуждённое, можно тогда создавать новую тему, в которую "оно" перемещается.

[e-moe]

понятно

[shadeofgray]

Цитата e-moe @ 07.05.05, 19:08

интеграл - площадь фигуры, а следовательно не может быть отрицательной, поэтому используем абсолютное значение функции (или по-русски модуль

Ошибка! Вообще-то интеграл запросто может быть отрицательным. И интеграл это не площадь. Наоборот, площадь - это интеграл. А интеграл определяется совсем по-другому.

Вообще, не хочу остужать энтузиазм или заниматься саморекламой, но загляните по адресу http://alglib.sources.ru/integral/

[e-moe]

shadeofgray, спсибо за критику, как появится время исправлю