Обещанное "ретуширование": символ ` воспринимайте как пробел.
Начинаю с ОТВЕТА (параметров этой формулы):
`u = 0,48892 68569 74369 ` А = 0,45416 39606 86749
v1 = 0,69088 05504 86348 `B1 = 0,21420 03609 26862
v2 = 0,93965 52580 96846 `B2 = 0,04273 12318 65773
wx = 0,91862 04410 56722 ``C = 0,14445 22232 60307
wy = 0,34487 20253 64404
Вычисления велись на 16-значном калькуляторе. Погрешности (± 1 или чуть больше) возможны лишь в последних знаках после запятой.
Первая строка означает, что при численном интегрировании по нашему квадрату функции f(x,y) следует сложить взятые с коэффициентом А величины f(u,0), f(0,u), f(-u,0) и f(0,-u).
Вторая и третья строки, - что надо продолжать складывать с коэффициентом В величины f(v,v), f(-v,v),
f(-v,-v) и f(v,-v). Здесь будет две такие группы.
Четвёртая и пятая строки, - что надо складывать с коэффициентом С восемь величин: f(wx,wy), f(wy,wx), f(-wx,wy), f(-wy,wx), f(-wx,-wy), f(-wy,-wx), f(wx,-wy) и f(wy,-wx).
Формула точна (если игнорировать погрешности в последних знаках) для всех многочленов до 9-й степени включительно.
Формула прошла тест-контроль при интегрировании одночленов:
1, `x^2, `x^4, `x^2 y^2, `x^6, ` x^4 y^2, `x^8, ` x^6 y^2, `x^4 y^4.
При этом отклонения от истинных теоретических значений интегралов наблюдались лишь в последнем (15-м) знаке после запятой. Симметричное расположение узлов кубатуры гарантирует аналогичные результаты при замене в этих одночленах x на y и наоборот. Кроме того, оно гарантирует равенство нулю численного значения интеграла от любого одночлена, содержащего хотя бы одно из переменных в нечётной степени, в частности, любого одночлена нечётной степени. Для выходящих за задуманный диапазон одночленов x^10, `x^8 y^2 и x^6 y^4 ошибки составили соответственно -0,74%, -0,17% и +3,7%.
- - - - - - - - - - - - - - - - - - - - -
А теперь хочу поделиться с неленивым (математически) читателем поразившими меня «чудесами». Для этого изложу (конспективно) вывод этой формулы.
Задача сводится к составлению и решению следующей системы девяти алгебраических уравнений с 9-ю неизвестными:
4 а1 `````+ 4 а2 ```` + 4 а3 ```` + 8 а4 ``````````````````````= 4
2 а1 x1^2 + 4 a2 x2^2 + 4 a3 x3^2 + 4 a4 (x4^2+y4^2) ``````````= 4/3
2 a1 x1^4 + 4 a2 x2^4 + 4 a3 x3^4 + 4 a4 (x4^4+y4^4) ````````` = 4/5
```````````4 a2 x2^4 + 4 a3 x3^4 + 8 a4 (x4^2 y4^2) ``````````= 4/9
2 a1 x1^6 + 4 a2 x2^6 + 4 a3 x3^6 + 4 a4 (x4^6+y4^6) ``````````= 4/7
```````````4 a2 x2^6 + 4 a3 x3^6 + 4 a4 (x4^4 y4^2+x4^2 y4^4) = 4/15
2 a1 x1^8 + 4 a2 x2^8 + 4 a3 x3^8 + 4 a4 (x4^8+y4^8) ``````````= 4/9
```````````4 a2 x2^8 + 4 a3 x3^8 + 4 a4 (x4^6 y4^2+x4^2 y4^6) = 4/21
```````````4 a2 x2^8 + 4 a3 x3^8 + 8 a4 (x4^4 y4^4) ``````````= 4/25
Избыток чётных коэффициентов связан здесь с самостраховкой: нет ли ошибок при составлении самой системы?
Здесь х1 – координата узла, расположенного на положительной части оси «х».
а1 – соответствующий этому узлу весовой коэффициент (используемые здесь
обозначения отличаются от обозначений в ОТВЕТЕ).
х2 и х3 – координаты узлов на главной диагонали: (х2,х2) и (х3,х3): ``` 0 < x2 < x3.
а2 и а3 - их весовые коэффициенты.
х4 и у4 – координаты узла (х4,у4); ```х4 > y4 > 0.
а4 – его весовой коэффициент.
Как видим, задача похожа на школьную (для старших классов). Однако сомневаюсь, что найдётся школьник, у которого хватит терпения довести её до конца…
Введём для сокращения писанины обозначения:
x1^2 = p ``````` 2 a1 = u
x2^2 = q ````````2 a2 = v
x3^2 = r (r > q) ``2 a3 = w ````````````` (0)
x4^2 = s ````````2 a4 = x
y4^2 = t (s > t)
Подставив это в систему и произведя естественные упрощающие вычитания уравнений, найдём, что из 8-го уравнения можно выразить величину
x = 8/(525st(s-t)^2). ```````````````````(1)
Если подставить это в остальные уравнения и изменить их порядок, получим систему из 8-ми уравнений с 8-ю неизвестными:
````2vq^4 + 2wr^4 `+ 32st/(525(s-t)^2) ``` = 4/25
````2vq^3 + 2wr^3 `+ 16(s+t)/(525(s-t)^2) ` = 4/15
````2vq^2 + 2wr^2 `+ 32/(525(s-t)^2) ``````= 4/9
up + 2vq ``+ 2wr ```+ 16(s+t)/(525st(s-t)^2) = 4/3
2u + 2v ```+ 2w ````+ 32/(525st(s-t)^2) `````= 4
up^2 ````````````` + 16/(525st) ``````````` = 16/45
up^3 ````````````` + 16(s+t)/(525st) ```````= 32/105
up^4 ````````````` + 16(s+t)^2/(525st) `````= 64/225
Последние три уравнения содержат 4 неизвестных. В результате не очень простого (но и не чрезмерно сложного!) анализа нужные нам комбинации можно выразить через величину p:
.
u = 256/(315p^2 (35p^2-60p+28))
s+t = (2/5) (15p – 14)/(7p – 6) ``````````````````````` (2)
st = (3/25) (35p^2 – 60p +28)/ (7p-6)^2 ```````````Равенства эти легко проверяются
(s-t)^2 = (32/25) (15p^2 – 30p + 14)/ (7p-6)^2 `````подстановкой в эти 3 уравнения.
Подставляя эти выражения в первые пять уравнений и перенося члены, содержащие величину p, в правые части, получим в результате весьма утомительных выкладок систему 5-ти уравнений с 5-ю неизвестными:
vq^4 + wr^4 = A
vq^3 + wr^3 = B `````````````````(3)
vq^2 + wr^2 = C
vq `` + wr ``= D
v ````+ w `` = E,```` где правые части приведены к общему знаменателю:
630p^2 (15p^2-30p+14) (35p^2-60p+28) ,
а числитель А = (693p^2 – 1440p + 655,2) (35p^2 – 60p + 28)p^2
- - - - - - - B = (945p^2 – 1956p + 924) (35p^2 – 60p + 28)p^2
- - - - - C = (1365p^2 – 2940p + 1420) (35p^2 – 60p + 28)p^2
- - - D = (91875p^5 – 368200p^4 + 546000p^3 – 363840p^2 + 96720p – 3584)p
- E = 361375p^6 –1428000p^5 +2091600p^4 –1360800p^3 +324240p^2 +15360p –7168
Кто-нибудь чувствует «родство» между этими 5-ю многочленами? Я во всяком случае – первоначально не чувствовал… К тому же числители D и E получились каждый в результате приведения к общему знаменателю и сложения трёх дробей с разными знаменателями. «Чудом» оказалось то, что оба они делятся на (35p^2 – 60p + 28). Как это было обнаружено? Случайно заметил, что 91875 и 361375 делятся на 35, а 3584 и 7168 – на 28. Возникло непреодолимое желание произвести соответствующие деления «в столбик». И всё разделилось!.. Итак, предыдущие строки можно заменить такими:
где правые части приведены к общему знаменателю:
630p^2(15p^2 – 30p + 14)
а числитель А = (693p^2 – 1440p + 655,2)p^2
- - - - - - - B = (945p^2 – 1956p + 924)p^2
- - - - - C = (1365p^2 – 2940p + 1420)p^2
- - - D = (2625p^2 – 6020p^2 + 3180p – 128)p
- E = 10325p^4 –23100p^3 +11900p^2 –256
ЗАМЕЧАНИЕ о проведении «утомительных» выкладок. При плохой организации работы их практически невозможно осуществить, так как шагов много, а вероятность ошибок на каждом шаге довольно велика. Причём в основном «дурацких» ошибок: перепутать цифру или знак (плюс-минус), пропустить символ, написать его два раза, выбрать число из «чужой» позиции… Бывают (реже) и более серьёзные ошибки, связанные со сбоем в понимании своей «теории».
Пример. Если вероятность ошибки на каждом шаге равна 0,01, то вероятность НЕ
ошибиться за 100 шагов равна (1-0,01)^100=1/e = 0,368, за 200 шагов 1/e^2 = 0,135, а за 300 шагов 1/e^3 = 0,05…
Подобно альпинистским восхождениям здесь не обойтись без страховки. Моя страховка состояла в следующем. Каждый ЭТАП вычислений состоит из серии элементарных шагов. Сделав один шаг, СРАЗУ делал его повторно и сравнивал результаты. Завершив этап, производил его численную проверку: придавал переменным конкретные числовые значения и сравнивал числовые значения преобразуемых выражений в начале этапа и в его конце. Иногда использовал более одного варианта придаваемых значений. Этот тест-контроль хотя и не является доказательством безошибочности выкладок, но значительно повышает уверенность и смелость продолжать работу. При этом объём труда возрастает всего примерно в 3 раза по сравнению с бесконтрольной работой, которую, как я сказал выше, практически невозможно завершить без ошибки.
ТЕОРЕТИЧЕСКОЕ ОТСТУПЛЕНИЕ.
В систему (3) неизвестные v и w входят линейно. Выразим v по «правилу Крамера» cперва из 1-го и 2-го, а затем из 2-го и 3-го уравнений. Получим два разных выражения, приравняв которые, получим после упрощений первую строку такой системы:
Сqr – B(q+r) + A = 0
Dqr – C(q+r) + B = 0 ````````````` (4)
Eqr – D(q+r) + C = 0
Вторая строка получается аналогично из 2-го и 3-го, а затем 3-го и 4-го уравнений.
Третья - - - - - - - - - - - - - - - - - - - - 3-го и 4-го, - - - - - 4-го и 5-го - - - - - - -
Система (4) выглядит как однородная система линейных уравнений относительно неизвестных: qr, (q+r) и z, ИМЕЮЩАЯ ЗАВЕДОМО ненулевое решение: (qr, (q+r), 1).
А это возможно лишь тогда, когда определитель такой системы равен нулю. Итак, НЕОБХОДИМЫМ условием разрешимости системы (4) является равенство:
C^3 + B^2 E + D^2 A – ACE – 2BCD = 0
Это – уравнение для нахождения величины p. Дальнейшие ещё более утомительные выкладки приводят к тому, что оно эквивалентно такому:
(155p^2-150p+27)/(p^2 (15p^2-30p+14)) = 0 ``````````````(5)
То есть задача свелась к решению квадратного уравнения!
В ходе этих выкладок было обнаружено (примерно аналогичным способом) ещё более удивительное «чудо»: делимость некоторого промежуточного многочлена 6-й степени на (15p^2-30p+14)^2 , что и привело к квадратному уравнению.
Решая (5) получим для p два значения:
p1 = (75-12√10)/155 = 0,23904 94714 708352…
p2 = (75+12√10)/155 = 0,72869 24640 130354…
Отсюда в частности видно, как получить параметры нашей кубатурной формулы с произвольным числом знаков: надо точнее вычислить √10 и все последующие квадратные корни (кроме этого будут использоваться лишь “четыре арифметических действия»).
Дальнейший план простой. Надо для p1 и p2 вычислить A, B, C, D и E. Затем решить систему (4). При этом обнаружится, что условие равенства нулю рассмотренного определителя не только необходимо, но и достаточно (окажется, что z не равно нулю!) для её разрешимости. Однако величину p2 придётся отвергнуть, т.к. при этом будет qr < 0. Затем, зная qr и q+r. найдём сами q и r (q < r). Потом из системы (3) по «правилу Крамера» найдём v и w.
Затем из соотношений (2) найдём u, st и s+t. Потом найдём s и t и, наконец, из соотношения (1) получим x. Осталось из соотношений (0) вычислить нужные узлы и весовые коэффициенты.
[свернуть]