Подынтегральная функция методом Симпсона

[Alessandro4]

Подынтегральная функция методом Симпсона

Помогите с задачей и объясните ее доступным для новичка языком:tehnari_ru_674:

Задача: вычисление определённых интегралов и табулирование первообразных функций

Все что мне известно:

Подынтегральная функция - 1/x*ln*x
Промежуток интеграла - [2;3]
Кол-во частей разбивки - 36
Вычисление значения первообразной h - 0,2
Точность вычисления значения первообразной - 0,001
Точное значение первообразной - 2,3026(ln ln x - ln ln 2)

 

[Vladimir_S]

Да легко, сейчас нарисую.

 

[Vladimir_S]

Значит, так.
Прежде всего — методичка составлена бестолково и ошибками. В частности, откуда взялся множитель 2,3026 — совершенно не представляю. А так — пожалуйста:

Const
 a=2.0;
 Eps=0.001;
 h=0.2;

Var
 x:Real;
 j:Byte;

Function F(p:Real):Real;
begin
 F:=1.0/(p*ln(p));
end;

Function Antiderivative(p:Real):Real;
begin
 Antiderivative:=Ln(Ln(p));
end;


Function Simpson(p:Real):Real;
var
 i,N:Integer;
 Sum1,Sum2,t:Real;
begin
 N:=2;
 t:=(p-a)/(2*N);
 Sum1:=(F(a)+F(p))/3*t;
 for i:=1 to N do
  Sum1:=Sum1+4.0*t/3.0*F(a+t*(2*i-1));
 for i:=2 to N do
  Sum1:=Sum1+2.0*t/3.0*F(a+t*(2*i-2));
 Sum2:=Sum1;
 Repeat
  Sum1:=Sum2;
  N:=N*2;
  t:=(p-a)/(2*N);
  Sum2:=(F(a)+F(p))/3*t;
  for i:=1 to N do
   Sum2:=Sum2+4.0*t/3.0*F(a+t*(2*i-1));
  for i:=2 to N do
   Sum2:=Sum2+2.0*t/3.0*F(a+t*(2*i-2));
 Until ABS(Sum2-Sum1)<Eps;
 Simpson:=Sum2;
end;

Begin
 WriteLn('   x     Simpson    Accurate value');
 for j:=1 to 5 do
  begin
   x:=a+h*j;
   WriteLn(x:5:1, '   ',Simpson(x):0:8, '   ',(Antiderivative(x)-Antiderivative(a)):0:8);
  end;
 ReadLn
End.

 

[Alessandro4]

Спасибо вам)) Я как дома буду проверю
А как вам мой вариант решения задачи, помог препод:

 

[Vladimir_S]

А как вам мой вариант решения задачи

Извините, но вижу лишь обрывок весьма, мягко говоря, странного кода. Симпсоном тут и не пахнет.

Я как дома буду проверю

А что проверять? Я ж выложил результат работы программы.

 

[Alessandro4]

Большое вам спасибо!
Я больше имел в виду что изучу ее, чтоб суметь грамотно объяснить преподавателю

 

[Vladimir_S]

Я больше имел в виду что изучу ее, чтоб суметь грамотно объяснить преподавателю

Удачи! И обращайтесь, если потребуются с моей стороны какие-то объяснения.

 

[Alessandro4]

Ладно, пытался я значит сам разобраться, но что то не очень, выходит:tehnari_ru_325:
Объясните, пожалуйста, ваши блоки и сам метод)

 

[Vladimir_S]

Объясните, пожалуйста, ваши блоки и сам метод)

Пожалуйста.
Интеграл берется аналитически, и первообразная (кстати, в методичке, как обычно бывает, это слово употреблено безграмотно, на самом деле первообразная, она же "неопределенный интеграл", это функция, дифференцирование которой даёт подынтегральное выражение, поэтому первообразная определяется с точностью до произвольной константы) есть
Ln(Ln(x))+С.
Отсюда определенный интеграл с переменным верхним пределом (который эти умники называют "первообразной") есть
Ln(Ln(x)) - Ln(Ln(a)).
Уж откуда эти ... взяли числовой коэффициент — не ведаю. Впрочем, догадываюсь: скорее всего, подынтегральная функция предполагалась не 1/(xLn(x)), как у них в таблице, а 1/(xLg(x)), т.е. с десятичным, а не натуральным логарифмом. Тогда действительно возникает коэффициент Ln(10). Идиоты!
Всё это считается в программе через функцию Antiderivative.

Теперь, собственно, о численном интегрировании.

1. Вспомним формулу Симпсона для вычисления интеграла на промежутке (a,b):
а) разбиваем промежуток на 2N частей (тут важна чётность).
б) вводим
h = (b - a)/(2N)
Xj = a + h*j, j = 0, 1, ... 2N
Yj = f(Xj), где f — подынтегральная функция.
тогда интеграл есть
(h/3)*(Y0 + 4*Y1 + 2*Y2 + 4*Y3 + ... + Y2N)

2. Тут опять "уважаемые" преподаватели бьют бедных студиозусов ключом (разводным или шведским) по голове. Видите ли, можно либо зафиксировать число разбиений интервала, тогда точность — уж как получится, либо задать точность и дробить интервал до тех пор, пока мы этой точности не достигнем. А то и другое сразу, как предписывают эти корифеи, НЕВОЗМОЖНО. Поэтому я выбрал второй вариант — вычисление интеграла с заданной точностью. Вычисление делается в функции Simpson.

3. В этой функции организован следующий цикл:
а) Вычисляется значение интеграла при 2N = 2. Обозначим это значение Sum1.
б) Удваиваем N и повторяем вычисление. Получаем величину Sum2.
в) Если |Sum1-Sum2|< ε, то на этом заканчиваем, в противном случае величине Sum1 присваиваем значение Sum2, снова вдвое дробим интервал и повторяем вычисление, в результате которого получаем новое значение Sum2, которое сравниваем с Sum1.
г) Делаем так до тех пор, пока не выполнится условие малости разницы результатов двух итераций.

Далее в программе сделан вывод таблицы результатов, где сопоставляется значение численного расчета при изменении верхнего предела интеграла с точным значением.

 

[Alessandro4]

 begin
1  N:=2;
2  t:=(p-a)/(2*N);
3  Sum1:=(F(a)+F(p))/3*t;
4  for i:=1 to N do
5   Sum1:=Sum1+4.0*t/3.0*F(a+t*(2*i-1));
6  for i:=2 to N do
7   Sum1:=Sum1+2.0*t/3.0*F(a+t*(2*i-2));
8  Sum2:=Sum1;
9  Repeat
10   Sum1:=Sum2;
11   N:=N*2;
12  t:=(p-a)/(2*N);
13   Sum2:=(F(a)+F(p))/3*t;
14  for i:=1 to N do
15   Sum2:=Sum2+4.0*t/3.0*F(a+t*(2*i-1));
16   for i:=2 to N do
17   Sum2:=Sum2+2.0*t/3.0*F(a+t*(2*i-2));
18 Until ABS(Sum2-Sum1)<Eps;
19  Simpson:=Sum2;
20 end;

Спасибо за разъяснение, но многое мне, конечно, не понятно
Кстати и препод сказал, что вы ошиблись хех, да да... а именно в строчке - 7, опираясь на формулу сказал, что "-2" не должно там быть и попросил меня понять как работает весь этот цикл > объяснить ему > сдать, ибо мое объяснение его не особо удовлетворило. Объясните, пожалуйста, все строчки, если вас не затруднит, ибо в методе симпсона и в самих интегралах я не очень, желательно разжевав) Еще раз огромное спасибо)

 

[Vladimir_S]

Хорошо.
Прежде всего — НИКАКИХ ОШИБОК В КОДЕ НЕТ! И идёт Ваш препод, извините, лесом...
Теперь давайте разбираться.
Как я писал выше, формула Симпсона есть
(t/3)*(Y0 + 4*Y1 + 2*Y2 + 4*Y3 + ... + Y2N)
Давайте для простоты предположим, что интервал (a, b) мы разбили на шесть отрезков (повторяю, число отрезков должно быть ЧЁТНЫМ). Таким образом, по оси X мы имеем 7 точек: X0, X1, X2, X3, X4, X5 и X6. Соответственно, по оси Y имеем 7 значений: Y0, Y1, Y2, Y3, Y4, Y5 и Y6, равных f(X0), f(X1), f(X2), f(X3), f(X4), f(X5) и f(X6), причем X0=a, Y0=f(X0)=f(a) и X6=b, Y6=f(X6)=f(b), где f — подынтегральная функция. Величина t есть (b-a)/6. Теперь давайте искать сумму.
1. Определяем вклад концевых точек. Он есть
t/3*(Y0+Y6)=t/3*(f(a)+f(b)).
2. Добавляем вклад от точек X1, X3 и X5. Соответствующие им значения Y1, Y3 и Y5 идут с коэффициентом 4. Поэтому добавляем к сумме величину
t/3*(4Y1+4Y3+4Y5)=t/3*(4f(X1)+4f(X3)+4f(X5))
3. И, наконец, добавляем вклад от точек X2 и X4. Соответствующие им значения Y2 и Y4 идут с коэффициентом 2. Поэтому добавляем к сумме величину
t/3*(2Y2+2Y4)=t/3*(2f(X2)+2f(X4)).
Всё!
Теперь построим алгоритм перебора индексов при нахождении указанных сумм. Заметим, что
X1=X0+t=a+t
X2=X0+2t=a+2t
X3=X0+3t=a+3t
X4=X0+4t=a+4t
X5=X0+5t=a+5t
X6=X0+6t=a+6t=b
Таким образом, для X1, X3 и X5 имеем
X(2i-1)=a+t*(2i-1), i меняется от 1 до 3 (слева в скобках — индекс).
Для X2 и X4 имеем
X(2i-2)=a+t*(2i-2), i меняется от 2 до 3.
Вот, собственно, всё это и проделано в строках от 1 до 7 и от 11 до 17. И, ещё раз повторяю, НИКАКИХ ОШИБОК ТАМ НЕТ!!!
Дальше там организован цикл мельчения разбивки до достижения требуемой точности.
ВСЁ!

P.S. Возможно, Ваш высокоучёный препод не заметил, что в строке 6 цикл берется в пределах от 2 до N. Можно, конечно, взять пределы от 1 до N-1, тогда "-2" в строке 7 нужно убрать. Но это — исключительно дело вкуса, на результат не повлияет. Советую Вам ткнуть этого ученого мужа носом в шестую строку, чтобы осознал значения поставленных пределов и не вякал про "ошибки". Хех, да да...

 

[Alessandro4]

Благодарю! Вроде потихоньку вникаю)
Еще вопросики.. В прошлом разъяснении вы сказали, что h=(b-a)/(2N), а в коде у вас, да и в моей методичке h уже известен - 0,2, насчет t, должно быть это коэффициент перед "скобкой" в формуле Симпсона, раз он равен (b-a)/6, а в самом коде вы присваиваете ему (p-a)/(2*N), почему? И я так понял переменную b вы не ввели и почему вы не учли кол-во частей разбивки - 36?

 

[Vladimir_S]

Отвечаю в порядке поступления.

В прошлом разъяснении вы сказали, что h=(b-a)/(2N), а в коде у вас, да и в моей методичке h уже известен - 0,2

Так, не путаем, как говорится, "божий дар с яичницей". В задаче речь идёт о вычислении интеграла С ПЕРЕМЕННЫМ ВЕРХНИМ ПРЕДЕЛОМ, т.е. о построении зависимости интеграла от х — верхнего предела. Так вот, h — это ШАГ ИЗМЕНЕНИЯ ВЕРХНЕГО ПРЕДЕЛА, к разбиению интервала при численном интегрировании для каждого значения х эта величина отношения НЕ ИМЕЕТ. Шаг разбиения интервала от а до х у меня обозначен буквой t. Действительно, обычно эта величина в справочниках и пособиях обозначается, как h, отсюда некоторая путаница. Извините.

насчет t, должно быть это коэффициент перед "скобкой" в формуле Симпсона, раз он равен (b-a)/6, а в самом коде вы присваиваете ему (p-a)/(2*N), почему?

(b-a)/6 — это только в поясняющем примере, не более того. В задаче же речь идёт о вычислении интеграла с ПЕРЕМЕННЫМ верхним пределом. Этот переменный предел является аргументом функции Simpson (обозначен p). При обращении к функции на его место подставляется фактический параметр х.

И я так понял переменную b вы не ввели

Не ввёл. За полной ненадобностью.

и почему вы не учли кол-во частей разбивки - 36?

Я уже подробно отвечал на этот вопрос, см. #9, пункт 2. За идиотизм авторов методички, извините, не отвечаю.