Путенихин Петр Васильевич
Численное интегрирование
ББК 22.12; 22.313; 22.62; 430ф
УДК 512; 51-7; 51-71; 53.023; 531-4; 531.51; 52-423
Путенихин, Петр Васильевич.
П90 Численное интегрирование / П.В.Путенихин. -
Саратов: ООО "АМИРИТ", 2022. - 144 с., цв. илл.
ISBN 978-5-00207-049-7
Оглавление
1. Базовые положения численного интегрирования
2. Притяжение точки к дуге
3. Притяжение точки к прямой линии
4. Притяжение двух линий
5. Притяжение двух плоских зеркал
6. Притяжение точки к шару
Притяжение двух шаров
7. Притяжение двух дисков
Точка в плоскости диска
Точка над плоскостью диска
8. Притяжение двух стержней
Соосные стержни
Параллельные стержни
9. График параметрической интегральной функции
Выводы
ПРИЛОЖЕНИЯ
П1. Excel - процедура
П2. Загадочное уравнение
Ссылки
1. Базовые положения численного интегрирования
Интегрирование, получившее название численного (числового), является способом вычисления постоянного интеграла, вычислить который традиционным математическим, аналитическим способом невозможно, поскольку он либо не является табличным, либо его приведение к табличному невозможно или чрезвычайно трудоёмко.
В математике вычисленный определённый интеграл всегда является определённым числом. Напротив, неопределённый интеграл, являющийся основой определённого интеграла, - это всегда некоторая функция, уравнение, переменная. Главным объектом интегралов являются подынтегральные функции, которые могут быть функциями любого числа переменных и параметров.
Вычисление определённого интеграла может производиться как аналитически, так и численными методами. Аналитическое вычисление - это двухэтапный процесс. Сначала интеграл вычисляется как неопределённый, после чего подставляются все его параметры и переменные.
Численное интегрирование - это прямое, непосредственное вычисление интеграла. Из этого следует, что перед вычислением определённого интеграла методом численного интегрирования все эти переменные и параметры должны получить однозначные числовые значения.
Поскольку любая подынтегральная функция - это связь, по меньшей мере, двух переменных, создаётся впечатление, что кратные интегралы описывают многомерные пространства - двухмерное и выше. Однако это не совсем верно. Конечно, набор переменных любого интеграла образует соответствующее многомерное пространство, мерность которого равна кратности интеграла, однако само вычисленное значение интеграла всегда просто число. Независимо от кратности интеграла его вычисленное значение является скалярной величиной. Скаляром в рассматриваемых случаях мы считаем любую численную величину, как безмерную, так и имеющую какую-либо размерность: метр, кубометр, килограмм на квадратный сантиметр и так далее. Размерность этого скаляра определяется размерностью пределов интеграла и назначается, выбирается тем, кто сформировал этот интеграл. В пространстве многомерного, кратного интеграла его вычисленное значение является точкой, положение которой в этом пространстве предопределено смыслом, заложенным в интеграл его разработчиком.
Например, определённый интеграл одной переменной, одинарный интеграл - это, скаляр, точка, определяемая переменными двухмерного пространства. Нет никаких веских аргументов, чтобы задать координаты этой точки на осях интегральной плоскости.
Подынтегральная функция одинарного интеграла геометрически на интегральной координатной плоскости, то есть, плоскости образованной осями координат интеграла, изображается в виде графика, некой кривой. В этом случае интеграл от этой функции геометрически представляет собой площадь криволинейной трапеции под этой кривой, ограниченной осью абсцисс и пределами интегрирования. Значение интеграла в этом случае является, с одной стороны, безразмерным числом, точкой, а с другой - двухмерным многообразием - площадью, имеющей размерность квадрата расстояния. В обоих представлениях нет возможности привязать их к конкретной точке двухмерной системы координат, координат подынтегральной функции.
Если рассматривать решение как безразмерную точку, то для её отображения можно к двухмерной плоскости добавить дополнительную координату - ось значений интеграла. Такой подход к решению интеграла создаёт плоскость, параллельную другим осям плоского пространства, а это уже не точка, а массив точек. Очевидно, у такого изображения, размещения точки в пространстве нет явного математического или физического смысла. Другой подход - это привязка к решению сложного массива координат: две координаты - это пределы интегрирования, а третья "координата" - это отрезок функции между этими пределами. Такое решение также выглядит весьма необычно.
Принцип вычисление одинарного интеграла как площади под графиком подынтегральной функцией можно расширить и на интегралы большей кратности. В частности, двойной интеграл можно рассматривать как объём тела, образованного двухмерной поверхностью подынтегральной функции и плоскостями пределов интегрирования. Задать координаты этого объёма одним числом, тем более, значением интеграла нельзя. Ещё более сложная конструкция - тройной интеграл. Можно представить его решение как отображение массы тела переменной плотности. Формальный механизм вычисления кратных интегралов имеет некоторое сходство, объединяющее значение интеграла как скаляра, и это скалярное решение с многомерными объектами, многомерным пространством интегральных переменных.
Для вычисления одинарного интеграла существуют, по меньшей мере, три принципиально различных способа. Первый способ - это аналитическое вычисление новой функции - неопределённого интеграла. Этот способ позволяет получить решение в точном аналитическом виде, что позволяет в дальнейшем осуществлять строгие аналитические преобразования и исследование интеграла, если только при его вычислении не использовались упрощающие замены, фактически подменяющие исходную функцию её подобием.
Второй способ - это графическое разбиение площади на рисунке на элементарные участки и их пересчёт.
Третий способ подобен такому же разбиению на участки, но аналитически, с использованием уравнения интегрируемой функции. Буквально, интеграл вычисляется как сумма конечного, хотя и большого количества слагаемых. Этот процесс и обозначается термином численное (числовое) интегрирование.
С помощью компьютера подынтегральная функция вычисляется для ряда значений переменной, а вычисленные значения сразу же суммируются. Умножение на дифференциал переменной, очевидно, можно произвести как при каждом вычислении подынтегральной функции, так и после суммирования, что в некоторых случаях уменьшит время выполнения компьютерной программы за счёт уменьшения числа умножений.
Этот способ позволяет получить либо график, либо таблицу числовых данных интеграла. Повторим, значение интеграла всегда единственное число, величина которого от интегральных переменных не зависит. Зависит она только от пределов интегрирования и величин параметров. Интегралы с неизменными пределами интегрирования назовём определёнными параметрическими интегралами. В зависимости от конкретных значений параметров такой интеграл имеет множество числовых решений, которые образуют собственное многомерное пространство. Иначе такой параметрический интеграл можно назвать интегральной функцией или параметрической интегральной функцией. Такая функция описывает собственное многомерное пространство, мерность которых определяется числом параметров. От многомерных пространств кратных интегралов, параметрические многомерные пространства отличаются физической трактовкой. Например, параметры могут не быть пространственными измерениями, что несколько диссонирует с представлениями о многомерности. Численное интегрирование всегда даёт либо конкретное число, либо таблицу чисел, которые могут быть представлены в виде графика, графической функции. Всё это позволяют графические возможности современных компьютерных приложений. Полученная в итоге некая линия, график, предельно точно соответствуют интегралу функции. Очевидно, её непосредственное аналитическое исследование невозможно. Конечно, отдельные участки вычисленного интеграла можно аппроксимировать какими-либо элементарными функциями, но фундаментальную суть подынтегральной функции и её интеграла такая аппроксимация, разумеется, не отражают.
Тем не менее, графическое, числовое (табличное) интегрирование даёт немало ценной информации. Если не предполагается дальнейшее фундаментальное исследование собственно интеграла, то такие графики вполне самодостаточны.
Рассмотрим подробнее, как выглядит такое числовое, табличное интегрирование на примере одинарного интеграла. Возьмём функцию, например, cos(x). Строим её график. Теперь вычисляем последовательно величины cos(x)dx и находим возрастающую сумму: cos1+ cos2+ cos3... для каждого текущего значения х. Пары значений (сумма; х) изображаем в виде графика. Результатом и является график - уравнение интеграла:
Рис.1.1. Числовой интеграл функции cos(x)
На рисунке представлены три графика. Первый - одиночная тонкая линия - это исходная, интегрируемая функция - cos(x). Толстая линия рядом - это полученный в реальном вычислении её числовой интеграл - sin(x). Рядом с толстой линией интеграла для сравнения показана аналитическая линии функции sin(x), которая немного смещена вправо, чтобы не сливаться с интегральной линией.
Подобный алгоритм позволяет производить интегрирование любой как аналитической, так и аппроксимированной, табличной функции. На рисунке рис.1.2 приведён пример трапецеидального интегрирования. Площадь под функцией, подынтегральной функцией f(x) разбивается на элементарные участки, имеющие форму трапеций. Аналитически вычисляется площадь каждой такой элементарной ячейки, которые затем суммируются. Интегрирование состоит из следующих шагов:
Находим первое значение функции.
Находим второе значение функции для прироста аргумента
Запоминаем его для следующего шага
Находим полу-сумму запомненного и нового значений
Умножаем на приращение аргумента (оно неизменно)
Получено значение площади первого интервала
В цикле находим остальные и суммируем их.
Вычисленный интеграл является определённым, на некотором интервале, поэтому каждая промежуточная сумма является точкой графика (таблицы) интеграла.
Рис.1.2. Площадь интегральной трапеции
При уменьшении приращения (x2-x1) площадь прямоугольной трапеции предельно приближается к площади косоугольной трапеции, одна из сторон которой является отрезком линии графика.
Механизм численного интегрирования был использован в работе [1] при исследовании решений уравнения теории относительности Эйнштейна. Рассматривалась задача о падении тела из состояния покоя на бесконечности на Чёрную дыру, имеющую заряд. Использовано решение Райсснера-Нордстрёма, преобразованное к интегралу следующего вида
Аналитическое интегрирование вряд ли возможно, поэтому для построения геодезической линии падающего тела использовано численное, табличное интегрирование. Для наглядности приведён график подынтегральной функции, имеющей следующий вид
Рис.1.3. График подынтегральной функции [1, с.255, рис.4.32]
Графический итог численного интегрирования приведён на рисунке рис.1.4.
Рис.1.4. Геодезическая частицы, падающей на Чёрную дыру Райсснера-Нордстрёма для внешнего наблюдателя и двух значений заряда Q [1, с.257, рис.4.33]
Анализ данного интегрального решения привёл к интересным и даже странным выводам, в частности:
"...можно заметить, что геодезическая падающей частицы выглядит не менее странно, чем график подынтегральной функции. Сам факт наличия трёх несвязанных отрезков практически невозможно согласовать с тем, что это, вообще-то, движение одной и той же частицы" [1, с.256].
Построение графика подынтегральной функции, строго говоря, не обязательно. Особой информации этот график не несёт, а для кратных интегралов его построение затруднено. Кроме того, такой график в некоторых случаях выглядит довольно странно. Например, некий условный график числа π.
Рис.1.5. Условный график числа π
Правильнее назвать этот график графиком разрядности числа π. Уравнение этого графика имеет вид
Это дискретный график, значение в каждой точке которого равно числу π, округлённому до d знаков после запятой. График выглядит весьма странно, поскольку число π - это константа, имеющая единственное значение. Странно привязывать график к единственному значению, однако это возможно.
Рассмотрим другие примеры использования численного интегрирования при решении некоторых задач в области гравитации. Отметим, что большинство из вычисленных далее интегралов являются параметрическими, в которых параметром является дискретность этих интегралов. Это означает, что собственно набор численных значений интеграла хотя и представлен в графическом или табличном виде, но на самом деле смысл их совпадает со смыслом "Графика числа π". Двумя значимыми величинами интеграла являются его предельные значения: предельное значение интеграла с максимальной дискретностью и значение этой дискретности. Наряду с ними, некоторые интегралы в данной работе имеют обычный графический смысл, схожий со смыслом графика подынтегральной функции рисунка рис.1.3.
2. Притяжение точки к дуге
Согласно закону всемирного тяготения Ньютона в классическом виде сила притяжения двух точек массой m равна
Здесь R - это расстояние между точками. Растянем одну из точек в дуговую линию с центром, совпадающим с другой точкой. Радиус дуги, следовательно, равен R, а полный её угол зададим равным некоторому значению 2α.
Рис.2.1. Притяжение точки к дуге
Результирующая сила притяжения точки к дуге равна
Вычисляем табличный интеграл
Найдём два значения силы притяжения масс. Сначала пусть угол α = 0. Это значит, что дуга сжата в точку, вследствие чего задача сводится к классическому закону всемирного тяготения
Здесь мы учли известное отношение: синус малого угла равен этому углу. Теперь рассмотрим другой пример. Пусть угол α = π/4. Получаем
Оказалось, что сила притяжения точки и дуги такой же массы составила примерно 90% от силы притяжения двух равных точек. Теперь для сравнения вычислим эту же силу по уравнению (2.2) численным интегрированием. Поскольку нас интересует разница значений, то вычисления произведём только для множителя с интегралом
Численное интегрирование на самом деле является простым суммированием конечного числа членов ряда. Ряд образуется значениями подынтегральной функции для последовательного ряда аргументов. Очевидно, значение этой суммы зависит от дискретизации, то есть, от количества дискрет аргумента в области его определения. Чем больше таких дискрет, тем сильнее значения суммы приближается к значению интеграла [3]. Нам неизвестно наилучшее количество дискрет, поэтому будет вычислять интеграл, меняя это количество. Выбираем предельно широкий диапазон разбиения: от 100 до 10'000.
На следующем рисунке приведён график изменения суммы в зависимости от количества дискрет. Там же проведена линия, соответствующая аналитическому значению интеграла. Линия горизонтальная, поскольку аналитическое значение интеграла соответствует разбиению его на бесконечное количество отрезков.
Рис.2.2. Численное интегрирование с разным числом дискрет
Замечаем странную картину. Значения интегральной числовой суммы при малых значениях разбиения скачкообразно меняются. Вероятно, это связано с погрешностями вычисления соответствующих дискретных значений подынтегральной функции. Действительно, VBA-процедура интегрирования построена по принципу, который назовём математическим. Диапазон изменения переменной интегрирования взят непосредственно в дробных значениях π/4, а шаг интегрирования получен делением этого диапазона на число дискрет, следовательно, является дробной величиной с большим числом знаков после запятой. При численном интегрировании в цикле производится многократное умножение на эту величину, затем суммирование результата. Собственно, использованная VBA-процедура содержит следующий фрагмент, в котором итогом вычисления является переменная k_duga:
step_ф = (Pi / 2) / n
k_duga = 0
For ф_duga = -Pi / 4 To Pi / 4 Step step_ф
dk_duga = Cos(ф_duga) * step_ф
k_duga = k_duga + dk_duga
Next
Очевидно, что погрешность от округления всех суммируемых дробей накапливается, причём в разных вариантах: с превышением и с уменьшением итога. Видимо, как результат и получена такая волнообразная кривая рис.2.2. На больших значениях дискретности расхождения нивелируются за счёт большого числа слагаемых. Избежать погрешности удобнее всего заменой дробной переменной интегрирования на целочисленную. Далее эту процедуру мы рассмотрим подробнее, здесь же приведём только результат её применения.
Рис.2.3. Численное интегрирование с разным числом дискрет на предельно большом интервале дискретизации
Для демонстрации того, что интегрирование привело к получению точного значения мы приводим контрольную линию, асимптоту, имя которой в легенде содержит её числовое значение - это значение аналитического решения. На всём интервале дискретностей, график выглядит гладким. Даже при большой дискретности скорость компьютерного вычисления интеграла достаточно высока, поэтому мы задали дискретность чрезвычайно большой. Практически любое значение дискретности выше 1000 даёт нам высокую точность интегрирования.
Разница силы притяжения между двумя точками и между точкой и дугой демонстрирует интересное явление. Казалось бы, массы взаимодействующих тел в обоих случаях одни и те же, а расстояние - как и между эквивалентными точками. Но силы - разные. Согласно закону тяготения, сила притяжения зависит от расстояния между телами. В случае с точками всё просто: расстояние между ними - это расстояние между их геометрическими центрами или центрами тяжести. Но если использовать тот же закон (2.1) к случаю точки и дуги, нам в это уравнение следовало бы подставить какое-то иное значение расстояния.
Здесь мы обозначили символом Rg радиус, условное расстояние между точкой и дугой, на котором сила их притяжения равна силе притяжения точек на расстоянии R. Этот радиус мы называем гравитационным. Преобразуем и находим
Мы нашли, что на всех дискретностях численного интегрирования и в аналитическом решении k < 1, следовательно, гравитационный радиус больше классического радиуса Ньютона
Буквально это можно трактовать как эквивалентное смещение центра тяжести дуги вдаль, за её пределы. Точка в центре дуги притягивается как бы к такой же точке, находящейся где-то снаружи, за пределами дуги. В частности, если дуга частично охватывает точку в центре, то гравитационный радиус может находиться сколь угодно далеко за пределами дуги.
В рассмотренной нами задаче притяжение точки, находящейся в центре дуги с углом π/2, составляет примерно 90% силы притяжения двух точек с такой же массой и расстоянии между ними, равном радиусу дуги.
3. Притяжение точки к прямой линии
Задача с дугой дала удобную возможность сравнить аналитическое вычисление силы притяжения с числовым интегрированием. При замене дуги на прямую линию такое сравнение становится невозможным.
Вытянем одну из точек в линию, длиной h. Оставшаяся точка находится над центром этой линии, на удалении d от неё. Линейная массовая плотность линии при этом составит
Для определения силы гравитационного притяжения объектов выберем на линии некоторый элементарный отрезок, точку dm на линии.
Рис.3.1. Гравитационное притяжение точки к линии
Отрезок имеет массу dm и притягивается к внешней точке m с силой dF
Проекция этой силы на вертикальную ось равна
Сила притяжения точки ко всей линии равна
Предполагаем, что сила притяжения точки к линии меньше силы притяжения точек, то есть, k < 1. Найдём этот уменьшающий коэффициент, приняв, h = 10, d = 0,001
Свести этот интеграл к табличному не удалось, поэтому вычисляем его методом числового интегрирования. Меняя значение дискретности числового интеграла n от некоторого минимального значения до осмысленно большого, находим соответствующие значения коэффициента.
Рис.3.2. Величина коэффициента уменьшения силы притяжения точки к линии в зависимости от дискретности численного интеграла
Как видим, график асимптотически приближается к величине
Решение оказалось предельно, убедительно точным. Смысл решения состоит в том, что точка притягивается к линии длиной h с силой, многократно меньшей, чем притяжение к такой же точке на таком же расстоянии. Если вновь использовать понятие гравитационного радиуса, то такая величина коэффициента означает, что точка притягивается к такой же точке, но находящейся на удалении в 200 раз дальше, чем прямая линия.
Другими словами, растягивание точки в линию привело к её эквивалентному гравитационному удалению в 200 раз.
4. Притяжение двух линий
Рассмотрим теперь гравитационное взаимодействие двух одинаковых линий. Для этого точки, которые изначально испытывали силу гравитационного притяжения (2.1), вытянем в линии, отрезки длиной h. Линейная массовая плотность линий при вытягивании составит
Выберем на линиях, которые обозначим A и B, элементарные отрезки, "точки" dx и dy
Рис.4.1. Гравитационное притяжение двух линий
Эти дифференциальные отрезки, точки имеют массу dm каждая и притягиваются друг к другу с силой dF
Проекция этой частной силы на вертикальную ось определяется углом наклона соединяющей линии и равна
Или окончательно, компактно
Полную силу притяжения находим интегрированием
Очевидно, что и в данном случае k < 1, то есть, сила притяжения линий друг к другу меньше силы притяжения двух m-точек в k раз. Найдём этот уменьшающий коэффициент
Для удобства разделим переменные
Интеграл такого вида ранее мы уже определили как не табличный и не сводимый к табличному, поэтому вновь его вычисление будем производить численным методом. Компьютерную форму интеграла удобнее представить с новыми переменными, дискретностью, изменяющейся в произвольно задаваемом диапазоне. Непрерывную, гладкую переменную x заменяем новой дискретной переменной u, ей эквивалентной. Если переменная x изменялась в интервале от 0 до h , то переменная u должна меняться в интервале от 0 до некоторого предельного дискретного значения n, определяемого дискретностью элементарных отрезков, то есть, их количеством в длине отрезка. В отличие от интегрального дифференциала, элементарные отрезки имеют конечное значение
Вводя новую переменную, мы устанавливаем между нею и заменяемой следующие соотношения
Фактически новая переменная обозначает количество элементарных отрезков в некоторой фиксированной длине исходной переменной. При изменении новой переменной в пределах количества элементарных отрезков n, исходная переменная изменяется в пределах полной длины отрезка. Проверка показывает корректность подобранной замены
Рассматривая связь между переменными как функциональную, мы нашли связь и между их дифференциалами. Подчеркнём: в компьютерных вычислениях дифференциалы являются эквивалентами конечных отрезков. Аналогично производим замену переменной y на переменную v.
Проверка показывает корректность и этой подобранной замены
Подставляем новые переменные в интеграл (4.1)
Преобразуем и сокращаем
Сделаем подстановку, упрощающую запись уравнения
Интегрирование показывает, что сила падает по мере увеличения дискретности интеграла. Однако при получении максимально точного результата возникла проблема: время интегрирования стремительно возрастало по мере увеличения дискретности n. Причём на первых значениях дискретности была обнаружена закономерность: увеличение в 2 раза дискретности приводило к такому же двукратному уменьшению прироста значения интеграла, что заметно в первых двух колонках трёх нижних строк таблицы, взятых в красную рамку. Фактически это означало отсутствие асимптоты, то есть, выглядело как бесконечное уменьшение интеграла вплоть до нуля. Понятно, что это плохой, неприемлемый результат.
Таб.4.1. Таблица значений числового интеграла (4.2)
Как выяснилось, эта закономерность была следствием недостаточно точного представления этих отношений в таблице. Поскольку интуитивно всё-таки предполагалось наличие асимптоты, предела значения интеграла, в таблицу и была добавлена колонка тенденции - колонка, в которой приведены более точные значения отношений значений интегралов двух смежных значений дискретности. То есть, отношение значения интеграла, вычисленного с некоторой дискретностью, к значению интеграла, вычисленного с удвоенной дискретностью. Эта колонка в таблице - крайняя правая. В результате асимптота, хотя и слабая, сразу же обнаружилась, поэтому было решено продолжить вычисления с более высокими значениями пределов, дискретности. По данным таблицы таб.4.1 построена диаграмма.
Рис.4.2. График зависимости коэффициента (синий) от дискретности интеграла и график тенденции (зелёный) его уменьшения
Графики являются функциями, зависимостями от дискретности n, представленной в логарифмическом виде log2n. Два основных графика на диаграмме имеют обозначения, имеющие в названии индекс "min". Следом за знаком равенства идёт минимальное значение этого графика на диаграмме. Основным, главным графиком является синий график коэффициента - "k_min=0,0199998". Его значения в таблице получены умножением на некоторый масштабный коэффициент значений интегралов. Масштабирование позволило избавиться на диаграмме от лишних нулей после запятой. Значения этой колонки и выведены на диаграмме. Заметим, что теоретическим значением, очевидно, должно быть 0,0200000..., поэтому на диаграмме приведена соответствующая этому графику асимптота со схожим названием "0,0200".
Верхняя зелёная линия, обозначенная в легенде как T_min = 1,00002, - это линия тенденции, то есть, напомним, отношения двух смежных значений интегралов. Тонкая линия, обозначенная в легенде как 1,0 - это асимптота, значение, к которому стремится в пределе линия тенденции. То, что отношение смежных значений интегралов в колонке тенденций и на диаграмме асимптотически стремится к единице, означает, что интеграл имеет конечное значение при достаточно большом количестве шагов интегрирования. В рассмотренном варианте, максимальное количество шагов интегрирования выбрано равным 409'600. Собственно метрическая дискретность линий, то есть, размер элементарных отрезков на каждой из них равен 10/409'600 ≈ 2,5·10-6 мм, то есть, почти в 106 раз меньше расстояния между линиями.
Вполне очевидным является вопрос: почему притяжение между двумя точками, растянутыми в линию меньше, чем между исходными точками? При этом уменьшение силы весьма существенно - в миллионы раз. Действительно, значение 0,02 интегралов и асимптоты в таблице и на диаграмме не учитывают масштаба. С учётом масштаба точное значение минимального значения коэффициента k_min и его асимптоты равно 2·10-5.
Однако ничего странного в таком уменьшении сил притяжения нет. Покажем это на простом примере. Рассмотрим две точки. Сила их взаимного притяжения равна (2.1):
Разделим эти две точки на четыре пары точек, то есть, каждую точку разделим на четыре одинаковые части. Сила притяжения каждых двух новых пар точек теперь будет в 16 раз меньше
Однако каждая из точек при этом притягивается и к смежным точкам из других пар. Непосредственно друг против друга точки размещены только в четырёх парах. Но каждая из точек испытывает притяжение ещё к трём смежным, лежащим в параллельной плоскости. Следовательно, суммарная сила притяжения будет состоять из 16 частных сил притяжения. Казалось бы, сила притяжения останется прежней
Однако это не так, поскольку лишь четыре силы из шестнадцати описываются уравнением (4.3). Это очевидно, поскольку сила притяжения частиц, не находящихся строго друг против друга, меньше силы, определяемой уравнением (2.1). Поскольку расстояние между этими "диагональными" частицами больше R, то и сила притяжения их будет меньше. Таких пар частиц - 12. Следовательно, общая сила притяжения равна
Величина в квадратных скобках меньше 16. Действительно, поскольку последняя дробь меньше единицы, то вся величина меньше 16.
Понятно, что это отклонение тем больше, чем дальше по горизонтали разнесены друг от друга все эти пары.
Для численного интегрирования всем параметрам мы присвоили конкретные числовые значения. Две точки, масса которых в вычислениях не учитывалась, находятся на расстоянии между ними d = 0,0001 мм. Растянув эти точки в отрезки длиной h = 10 мм, мы нашли, что сила гравитационного притяжения между линиями оказалась меньше силы притяжения между точками c учётом масштаба в k ~ 0,02·10-6 = 2·10-8 раз. В терминах гравитационного радиуса это уменьшение силы эквивалентно разнесению исходных точек на расстояние, примерно в 104 раз больше исходного.
5. Притяжение двух плоских зеркал
Рассмотрим теперь силу гравитационного притяжения между двумя площадками, зеркалами, образованными "расплющиванием" двух точек. Сравним эту силу с силой гравитационного притяжения (2.1), которая возникает между исходными точками, масса которых равна массе зеркал. Зеркала представляют собой квадратные площадки со стороной h, бесконечно малой толщины и массой m. Закон Ньютона в классическом виде (2.1) в этом случае неприменим. Для определения силы притяжения такой распределённой массы разобьём каждое зеркало на n-частей, площадь каждой при этом зададим равной k2. Очевидно, n = S/k2 = h2/k2. Между массами, сосредоточенными в точку, возникает сила гравитационного притяжения (2.1)
Если же массы распределены, то есть, точки расплющены, раскатаны в площадки, то суммарная сила будет определяться взаимодействием каждой элементарной площадки одного зеркала с каждой элементарной площадкой другого зеркала
Рис.5.1. Сила гравитационного притяжения элементарных площадок на двух зеркалах
Сила гравитационного взаимодействия двух таких элементарных площадок dA и dB равна
Справа в уравнении (5.1) каждый дифференциал массы - это масса соответствующей элементарной площадки dA или dB. Сами дифференциалы масс dm и их произведение, соответственно, равны
Уравнение (5.1) с их учётом приобретает вид
Радиус RAB для каждой, любой пары площадок равен
Подставляем значение Rxy, которое определяем согласно рисунку
Сделаем запись более компактной. Для этого заменим обозначения координат в виде названий точек с индексами названиями самих координат, то есть, их индексами. Получаем
Однако эта сила направлена вдоль радиуса RAB, а нас интересует сила, направленная вдоль интервала d, то есть, сила, притягивающая зеркала друг к другу. Из геометрических соотношений видим, что эта сила вдоль d во столько же раз меньше силы вдоль RAB, во сколько d меньше радиуса RAB
Следовательно, искомая сила равна
Учитывая константы, интегрированием находим результирующую силу гравитационного притяжения между зеркалами
Раскрываем, разворачиваем интеграл по каждой переменной
Интеграл, очевидно, табличным не является. Поэтому его значение будем искать с помощью числового интегрирования. Повторим, в этом механизме вычисления интегрирование заменяется суммированием конечного числа слагаемых, что легко осуществить с помощью компьютера. В нашем случае эта сумма выглядит так
Поскольку эквивалент дифференциала Δx содержится в каждом слагаемом, мы вынесли его за знак суммы. Смысл интеграла (5.2) довольно прост. Внутренний интеграл по первой переменной x вычисляет суммарную силу гравитационного притяжения некоторой точки на верхнем зеркале к каждой точке на оси x. При первом проходе - это точка левого ближнего угла верхнего зеркала.
Второй интеграл по переменной y вычисляет последовательно такую же силу притяжения к этой же точке следующего ряда точек оси x, имеющих новое, одно и то же значение координаты y. Интеграл вычисляется многократно, каждый раз обращаясь к вычисленному внутреннему интегралу предыдущего цикла. Завершение цикла этого интеграла по переменной y даёт суммарную силу притяжения всех точек нижнего зеркала к некоторой, всё той же точке на верхнем зеркале.
Следующий интеграл, по переменной v при первом проходе производит вычисление суммарной силы притяжения ко всем точкам нижнего зеркала всех точек одноименной оси, ближнего края верхнего зеркала. Соответственно, последний, внешний интеграл по u добавляет к общей сумме силу притяжения всех v-линий верхнего зеркала. Заметим, что порядок интегрирования, последовательность переменных значения не имеет, приведённый использован просто для демонстрации механизма численного интегрирования.
Такое интегрирование предполагает некоторое конечное значение дифференциала и, следовательно, их конечное число. Предел интегрирования, равный h, предполагалось разбить на n = 1000 дифференциальных отрезков. Следовательно, общее число вычислений подынтегральной функции равно n4 = 1012.
Нас интересует, насколько эта сила отличается от силы по закону Ньютона (2.1). Найдём их отношение. Присвоим коэффициенту индекс, совпадающий с пределом интегрирования
Если это отношение меньше единицы, то сила притяжения плоских зеркал меньше, чем притяжение точек такой же массы. Удобство уравнения (5.3) заключается в том, что оно содержит только геометрические параметры системы - отсутствует гравитационная постоянная и массы. При необходимости значение силы мы можем в дальнейшем вычислить с помощью (2.1) и этого коэффициента (5.3).
То же самое в развёрнутом виде
Хотя определённый интеграл (5.3) выглядит довольно красиво, его вычисление, численное интегрирование связано с серьёзной проблемой. В проведённых компьютерных вычислениях мы заметили, что в интегральной сумме необходимо использовать довольно высокую дискретизацию, доходившую порой до миллионов. Собственно перебор дискретностей интеграла направлен на определение такой дискретности, при которой её дальнейшее увеличение не ведёт к изменению значения интеграла, что можно назвать схождением интеграла. Если предположить, что интеграл (5.3) сойдётся хотя бы при дискретности 100, то общее число вычислений подынтегральной функции составит 108 раз. Это не просто счётчик от 1 до ста миллионов - это вычисление довольно замысловатой функции сто миллионов раз. Предварительное интегрирование показало, что даже дискретность 400 не ведёт к схождению интеграла. А это, следует отметить, уже почти в 256 раз больше, то есть, двадцать пять миллиардов вычислений функции.
Интегрирование с дискретностью 400 на использованном, далеко не самом слабом компьютере в VBA Excel 2010 заняло время больше суток. Напомним: дискретностью исследуемых определённых интегралов мы называем величину
Таким образом, одной из серьёзнейших, главных проблем числового интегрирования является недостаточное быстродействие компьютера. В отношении кратных интегралов, выше двойных это обстоятельство является решающим фактором. Например, вычисление интеграла (5.3) на интервале от 0 до 100 шагов с дискретностью n = 800 длилось почти 10 часов. Это значит, что полное решение интеграла с этой дискретностью займёт более 80 часов. При этом нет ясности, даст ли эта дискретность предельное, минимальное значение коэффициента kh. Поэтому необходимо произвести вычисление интеграла с большей дискретностью, например, n = 1'600. Время интегрирования возрастёт, видимо, в 24 = 16 раз. Если разбить интеграл на интервалы по 100 шагов, то их будет 16, а время вычисления одного интервала составит не 10 часов, а 160 часов. Полное время интегрирования составит 2'560 часов. Понятно, что для наших исследований это неприемлемое время. Поскольку в нашем распоряжении нет суперкомпьютера, для ускорения вычислений внесём в уравнение интеграла (5.3) дополнительные изменения. Сначала заменим переменные и пределы интегрирования, а в подынтегральной функции заменим некоторые устойчивые группы параметров равными им условными переменными. Замену переменных проведём по примеру задачи о притяжении двух линий. Используем приведённое выше соотношение
Новые переменные обозначим теми же именами, но со штрихом. Связь всех переменных, согласно соотношению, имеет такой же вид, как, в частности, для переменной х
Прежняя и новая переменные должны изменяться в сопоставимых пределах
Фактически новая переменная обозначает количество элементарных отрезков на некоторой фиксированной длине исходной переменной. При изменении новой переменной в пределах количества элементарных отрезков n, исходная переменная изменяется в пределах полной длины отрезка h. Проверка показывает корректность подобранной замены
По аналогии находим соотношения для трёх других переменных
Подставляем новые переменные в уравнение (5.3), меняем пределы интегрирования и сразу отбрасываем штрихи. В связи с заменой пределов интегрирования также изменим на n и индекс у коэффициента
Сделаем ещё одно очевидное преобразование: вынесем за скобки дроби в знаменателе подынтегральной функции
Обозначим дробь d/h символом p
Видим, что значение коэффициента kn не зависит ни от площади зеркал, ни от расстояния между ними, не зависит от их абсолютных значений. Зависит коэффициент только от отношения этих двух величин. Далее для ускорения вычислений явочным порядком в компьютерной процедуре выносим за пределы внутренней подынтегральной функции по x вычисление упомянутой устойчивой группы параметров, не зависящей от x, которую условно назовём групповым параметром q.
Вычисление группового параметра вне внутреннего интеграла, цикла по x, очевидно, уменьшит время выполнения цикла, ускорит вычисление интеграла. Далее во всех интегралах нижние пределы мы меняем с нуля на единицу, либо увеличиваем их на единицу.
Связано это с особенностями использованной в процедуре численного интегрирования конструкции For...Next. Процедура анализирует все значения пределов, то есть, производит вычисления как с нижним, так и с верхним пределом. В этом случае получается, что элементарных отрезков по осям на один больше. Интеграл на интервале от 0 до 1, кажущийся единичным, содержал бы не один шаг, как ожидается, а два. Единичным интервал 0...1 является в математическом интеграле. Добавим, что математический интеграл - это сумма бесконечного числа значений. Разумеется, можно уменьшить на единицу верхний предел. Но это менее удобно, поскольку во всех выражениях придётся вместо верхнего предела n вписывать более длинное выражение - (n - 1). При этом в строковых позициях по-прежнему придётся указывать n, что может ввести в заблуждение.
Как предполагалось и определено в процессе пробных вычислений интеграла, время его вычисления при больших дискретностях неприемлемо велико. Поэтому было решено производить вычисления интервалами. Интервальное вычисление интеграла означает его разбиение на несколько интегралов с последовательными пределами, после чего эти части суммируются. Такие вычисления можно производить как в разные периоды времени, так и на разных компьютерах. Собственно численный интеграл - это сумма некоторого конечного числа вычисленных значений подынтегральной функции. Например, если мы произведём численное интегрирование не до верхнего предела, а до его половины, то мы получим сумму первой, начальной половины из общего числа слагаемых. Если далее произвести суммирование, интегрирование от этого половинного предела, взяв его в качестве нижнего, до прежнего верхнего, то мы получим вторую половину интеграла.
Строго говоря, исходя из представлений о смысле интеграла, очевидно, что сумма последовательных, без разрывов любого числа интервальных интегралов равна интегралу на полном интервале. В общем виде интеграл для каждого интервала имеет вид
Здесь Δ - это величина интервала, а s - номер этого интервала. Для интервальных вычислений в форме интеграла (5.4) заменены пределы только первого интеграла, то есть, интервалы подбираются только по последней, внешней переменной. Проверим это равенство явным образом для численного интегрирования, разбив полный интервал на четыре равные части, то есть, задаём smax = 4.
Поскольку время интегрирования на больших дискретностях чрезмерно велико, в целях тестирования рассмотрим минимальные дискретности: 50, 100, 200, 400 и 800. Интервалы, как видим, выбраны с двукратным увеличением. Для этих вычислений уравнение интеграла и параметры имеют следующий вид и значения
Подчеркнём: значение интеграла (5.5) - это часть полного значения интеграла, равная 1/4 его полного значения. Иначе говоря, для получения полного значения интеграла следует просуммировать четыре интеграла, вычисленные на этих четырёх интервалах. Вычисляем интеграл (5.5) для указанных значений дискретности и сводим результаты в таблицу таб.5.1.
Таб.5.1. Значение интеграла (5.5) на полном интервале равно сумме его частных интервалов
Сразу же по данным таблицы замечаем, что интеграл (5.5) обладает интересным и важным свойством: все его интервальные значения равны друг другу. Для наглядности значения из таблицы приведём в виде диаграммы рис.5.2
Рис.5.2. Интеграл на полном интервале равен сумме интервальных интегралов
На диаграмме синяя линия - сплошная, хотя выглядит штриховой. Это график непрерывных интегралов, на полном диапазоне дискретностей. Красная линия - штриховая, это такой же график, но построенный по суммам четырёх интервальных интегралов. Линия выбрана штриховой, поскольку она полностью совпадает с синей и, будучи сплошной, могла бы её закрыть. График и таблица демонстрируют высказанное утверждение о равенстве интеграла и суммы его интервалов. Зелёный график, тенденция - это отношение двух смежных значений интегралов.
К месту отметим обстоятельства, описанные выше, в первом разделе работы. На диаграмме представлены три графика (не считая асимптот). Но в смысле понятий численного интегрирования, все эти графики - эквивалентны трём конкретным, единичным числам. Графики являются механизмом подбора, выборки этих трёх чисел: значения интеграла (5.5), контрольного параметра - тенденции, являющегося индикатором того, что подбор ведёт к конечному, предельному значению интеграла. Третий график служит демонстрацией совпадения значений интеграла на полном интервале с суммой его интервальных значений.
Основной точкой такого сравнения являются их величины на предельном, максимальном значении дискретностей. Данный тест, сравнение мы провели для того, чтобы убедиться в возможности определения численным интегрированием значения полного интеграла по его интервальному значению. В таблице таб.5.1 и на диаграмме рис.5.2 отчётливо видно, что все вычисления можно было произвести лишь на одном коротком интервале, в данном случае, на четверти полного интервала, и получить полное значение простым умножением полученного результата. Понятно, что это сокращает время интегрирования, в данном случае - в 4 раза.
Однако и четверть интервала на больших дискретностях требует длительных вычислений. Закономерно сделать ещё одну проверку на возможность ускорения вычислений. Для этого в том же уравнении (5.5) задаём минимально возможный шаг, интервал, равный единице, что в нашем случае означает равенство пределов. Величина smax в этом случае приобретает смысл периодичности интервалов. Уравнение (5.5) для единичного интервала приобретает следующий вид
Мы разбили полный интервал на 20 частей, smax = 20. На границах интеграла решение имеет несколько спорный вид, поэтому просто "отодвигаем" пределы интегрирования от границ на 1 - 2 шага и добавляем нулевую, 21-ую точку. При такой замене для получения полного значения интеграла, умножить на дискретность следует полученное в (5.6) усреднённое значение интервалов, поскольку уравнение даёт значение одного шаг из полного числа дискретностей. Использование среднего значения вычисленных интервалов необходимо, поскольку они, строго говоря, не равны друг другу. Итоги интегрирования представлены в следующей таблице
Таб.5.2. Единичные интервалы интеграла (5.6) с высокой точностью равны друг другу
Каждая из колонок - это значения интеграла, вычисленного с пределами, заданными номером s из первой колонки. И вновь обращаем внимание: значения интегралов в пределах каждой колонки практически тождественные. Расхождение значений просматривается далее, чем в 7 - 8 знаке после запятой. Это весьма интересная особенность интеграла.
В таблице значение каждого из интегралов соответствует верхнему пределу и определяется его дискретностью. Как показано в уравнении (5.6), значение верхнего предела равно произведению номера строки s на шаг интервалов Δ.
Верхний и нижний пределы мы выбрали равными, что означает длину, интервал интегрирования, равный единице. Например, в строке s = 3 в колонке s200 приведено значение интеграла с дискретностью 200 и пределами, равными 30.
Соответственно, значение интеграла в строке 8 и колонке s800 получено на дискретности 800 и с пределами
Напомним, что главной целью вычислений является определение такой дискретности интегралов, при которой её увеличение не ведёт к изменению значения интеграла. Это значит, что отношение значений в соседних столбцах таблицы таб.5.2 должно быть равно единице. Однако сравнение показывает, что для приведённых колонок это не так. Реальное отношение полных значений интеграла в таблице представляет последовательность
Таб.5.3. При увеличении дискретности отношение двух смежных значений интеграла (5.6) стремится к единице
Поскольку каждая строка соответствует единичному интервалу числового интегрирования, полное значение интеграла равно произведению любой из строк на дискретность. Понятно, это приблизительное значение, что следует просто из различия интервальных значений. В колонке Тенденция приведено отношение значения интеграла в строке выше, с меньшей дискретностью, к интегралу в текущей строке, дискретность которого выше. Просматривается слабая тенденция к уменьшению отношения и, видимо, его асимптотическое стремление к единице, хотя рассмотренные дискретности от этой асимптоты явно далеки. С другой стороны, уже на этом этапе время интегрирования на дискретности 1600 даже для единичного интервала оказалось весьма значительным, почти час. То есть, для следующей дискретности в 3200 время вычисления такого же единичного интервала составит почти 8 часов. Среднее время вычисления каждой строки меняется между колонками в геометрической прогрессии. Для каждой следующей дискретности время возрастает в 8 раз. Эта величина, например, для дискретностей 400 и 800 определяется как
Здесь параметры n с индексами - дискретность, верхние пределы трёх внутренних интегралов в (5.6) для двух смежных колонок. Подчеркнём, что в этом уравнении для N сравниваются два "соседних" интеграла с разными дискретностями. Для наглядности приводим таблицу таб.5.2 в графическом виде. Сначала в полном виде с логарифмической шкалой значений интегралов
Рис.5.3. Диаграммы k(n) для дискретностей n = 100...1600 и единичных интервалов s% от дискретности. Шкала k(n) логарифмическая
Обращаем внимание: и в таблицах и на диаграмме значения интегралов уменьшаются при увеличении дискретности, то есть, ожидается некоторое минимальное, предельное значение интеграла. Значения в строках колонок различаются незначительно, поэтому графики на рис.5.3 получились чрезмерно, демонстративно условными, можно сказать, декоративными. Следует признать, что такая равномерность выглядит несколько странно. Для демонстрации того, что интервальные интегралы в колонках всё-таки не равны друг другу, воспользуемся следующим приёмом отображения. Отбрасываем слева у всех значений в колонках равные знаки, оставляя только различающиеся. Например, у значений в колонке s200 совпадают первые 8 знаков после запятой. Значение в строке 1, следовательно, будет переписано следующим образом
Этот приём можно назвать отсечкой, удалением постоянной составляющей. В таком варианте изменения, различия будут видны существенно лучше. Поскольку графики новых дробей могут оказаться на неопределённых расстояниях друг от друга, добавляем вручную к ним некоторое положительное смещение, чтобы на диаграмме максимально их сблизить. Произведя замены во всей таблице, строим семейство графиков, зависимостей интегралов от их номера s в процентах.
Рис.5.4. Диаграммы k(s%) без постоянной составляющей
Теперь неравномерность, изменчивость интервальных интегралов в колонках видна предельно отчётливо. Заметим, что интегралы этих графиков равны значениям соответствующих полных интегралов. Как мы и ожидали, таблица таб.5.2 и рис.5.3 демонстрируют интересную картину. Независимо от дискретности, при выбранном равенстве интервалов интегрирования все интервальные значения интеграла оказались равны друг другу с очень высокой точностью. Кроме того, произведение каждого интервального значения на дискретность с высокой точностью равна значению интеграла, вычисленного на полном интервале переменной - от 1 до дискретности. Например, для строки s10 таблицы таб.5.2, соответствующей 50%, и строки 50% таблицы таб.5.1 расхождение имеют весьма малые значения, много меньше 1%. Подобное равенство значений наблюдается и для других уравнений, значения которым сведены в таблицу таб.5.4
Таб.5.4. Сравнение значений величин разных интегралов
Каждая дискретность представлена двумя строками: абсолютными значениями интегралов и их отклонениями в колонках от базового значения - первой колонки. Первая колонка - это значения интегралов на полном интервале дискретности, вычисленные в уравнении (5.6). Вторая колонка - это значение интеграла, равное сумме интервалов в 25%, вычисленное также в уравнении (5.6). Третья колонка - это значения интегралов (5.6), равное усреднённому значению единичных интервалов, умноженное на соответствующую дискретность. И, что особенно интересно, четвёртая колонка - это значения единичных интервалов интеграла (5.7), умноженное на дискретность. Интеграл (5.7) является тройным интегралам, что кардинально отличает его от четверного интеграла (5.6). Переход от четверного (5.6) к тройному интегралу (5.7) стал возможен благодаря обнаруженному свойству равенства интервалов (5.6).
Следует отметить, что обнаруженная особенность интервального равенства оказалась, вообще-то, неожиданной, но весьма интересной и, вместе с тем, чрезвычайно полезной. Действительно, при её использовании нет необходимости вычислять интеграл на полном интервале, от 1 до величины дискретности. Достаточно вычислить его на интервале, равном, например, единице, после чего просто умножить полученное значение на величину дискретности. Поэтому эту обнаруженную особенность (5.6) мы будем использовать для вычисления нашего четырёхкратного (четверного) интеграла (5.4) на больших значениях дискретности. Мы исходим из разумного предположения, что эта особенность будет иметь место и на других, на больших дискретностях и интервалах. Для вычисления интеграла на таких дискретностях это, пожалуй, единственное доступное решение.
Как мы отметили, полное значение интеграла, согласно выявленной особенности, равно произведению интеграла для одного шага на число шагов. Внесём в уравнение (5.4) соответствующие изменения
Напомним, что нижний предел здесь и далее изменён с нуля на единицу, что связано, как мы выше отметили, с особенностями используемого в компьютерной программе цикла For...Next. Ещё раз подчеркнём: это уравнение интеграла не математическое, в котором итог на нулевом интервале будет равен нулю, а компьютерное. Интервал 1 - 1 в данном случае означает, что вычисление подынтегральной функции, внутреннего тройного интеграла будет произведено только один раз, поэтому знак внешнего интеграла можно опустить. Четырёхкратный (четверной) интеграл, таким образом, сводится к интегралу тройному
Значение интеграла будет, очевидно, приблизительным, неточным, но мы предполагаем, что погрешность будет в приемлемых границах. Вместе с тем, следует отметить, что причину такого выравнивающего поведения интервальных интегралов объяснить весьма непросто. И здесь возникает естественный вопрос. Если сокращением интервалов удалось исключить внешний интеграл, может быть и три остальные также можно свести к единичным вычислениям?
Первый раз мы понизили кратность интеграла, уменьшив пределы внешнего интеграла до единицы и умножив полученное выражение на дискретность остальных интегралов. Уравнение (5.7) мы получили именно таким способом. Заполнив таблицу таб.5.2 значениями интервальных интегралов, мы обнаружили, что оставшиеся тройные интегралы на единичном интервале также кратно меньше полного интеграла. Если установить внешнему, третьему интегралу верхний предел в единицу, не окажется ли он вновь ровно в n-раз меньше полного интеграла?
Более того, малые погрешности, отклонения при уменьшении интервалов интегрирования таб.5.4 подсказывают ещё одну гипотетическую возможность ускорения вычислений и даже аналитического предсказания значения асимптоты графика kn понижением кратности всех интегралов. Однако оказалось, что такой переход ведёт к неприемлемым результатам. Для проверки этого обратимся к уравнению (5.7), развернув обратно подстановочные и групповые параметры.
Изменяем интервал интегрирования внешнего интеграла на предельно короткий, задав пределы от 1 до 1. После такой подстановки интеграл (5.7) приобретёт вид
Алгоритм численного интегрирования предполагает, что теперь внешний интеграл вычисляется только один раз, то есть, его значение в точности равно подынтегральной функции
Из простых логических соображений, используем эти же рассуждения и делаем ещё два подобных шага, преобразования
И последнее преобразование
Однако это явно неверный результат: предыдущие вычисления свидетельствовали, что коэффициент существенно меньше единицы. Поскольку такое аналитическое "понижение" уровня интегралов неприемлемо, далее будем вычислять интеграл (5.7) непосредственно, как есть, постепенно повышая дискретность, верхние пределы.
Поскольку интервалы 100 и 200 явно не дают приемлемых значений интегралов, далее рассматривать их будем лишь при особой необходимости. Итоги вычисления интегралов по уравнению (5.7) для широкого диапазона дискретностей представлены в таблице таб.5.5
Таб.5.5. Вычисленные значения интеграла (5.7) уменьшаются, а его тенденция асимптотически "пилообразно" стремится к единице
Первая колонка n - это дискретность интеграла, величины верхних пределов. Смысл двух других колонок полностью определяется их названиями. Значения в колонке Интеграл приведены в масштабе, чтобы избавиться от экспоненциального представления. Для наглядности на следующем рисунке приведены график значений интеграла, коэффициента kn и график его тенденции
Рис.5.5. Диаграмма к таблице таб.5.5
Левая вертикальная ось диаграммы относится к графикам интеграла, коэффициента kn и асимптоты 0,098. Правая вертикальная ось диаграммы относится к графикам тенденции и асимптоты 1,000. Специфический "пилообразный" вид графика тенденции вызван тем, что графики приведены с неравномерной шкалой дискретностей, горизонтальной осью диаграммы. Если использовать равномерную последовательность, ось, то график станет гладким, но общее время вычислений за счёт дополнительных промежуточных вычислений существенно возрастёт. Следует отметить, что в таком варианте с неравномерной шкалой информативность графика довольно плохая, его интерпретация затруднена. Если оставить только равномерно возрастающие значения, выделенные в таблице таб.5.5 цветом, то форма графиков становится достаточно гладкой
Рис.5.6. Сглаженные диаграммы к таб.5.5
Мы обоснованно предполагаем, что, интегральный график коэффициента kn, колонка Интеграл обязан асимптотически стремиться, и по наблюдениям он действительно стремится к какой-то предельной величине. Как следствие, график отношения двух смежных значений интеграла, который мы назвали графиком тенденции, также должен асимптотически стремиться к предельному значению - к единице. График тенденции показывает, что при увеличении дискретности, изменение коэффициента, интеграла обязано стремиться к нулю: два смежных его значения должны быть равны. Однако на этапе дискретностей выше 3600 необходимое время вычислений неприемлемо возросло, и достичь этого соотношения не удалось.
Следует вновь искать другие подходы. Выше мы обнаружили полезную особенность: значение интеграла можно вычислить, умножив его значение на коротком интервале на число этих интервалов в полном интервале. Очевиден вопрос: а можно ли оставшийся тройной интеграл вычислить подобным способом, не используя единичные интервалы во внутренних интегралах?
Произведём пробные вычисления интеграла (5.7), с пределами внешнего интеграла, которым теперь уже является третий, разбитыми на четыре равные части, по 25% от общей длины. Уравнение (5.7) приобретает следующий вид, где параметр s = 1...4
Результаты вычислений приведены в таблице таб.5.6
Таб.5.6. Интервальные значения интеграла (5.8) равны друг другу независимо от дискретности
Видим, что в таблице значения в каждой из колонок s25...s100, умноженные на 4, равны значению колонки k(n). Сравнение приведено для наглядности в форме погрешностей, которые представлены в колонках D25...D100. Например, в колонке D25 приведено отношение колонки kn к колонке s25, за вычетом единицы. Реальное отношение имело бы вид 1,000000009, что при любом представлении числа привело бы к расширению таблицы.
В использованном варианте мы просто "гасим" значащую единицу, а оставляем только дробную часть в экспоненциальной форме. Такое отношение даёт нам компактную погрешность, расхождение значений
Такое отношение в экспоненциальном виде максимально информативно и компактно, что, несомненно, более удобно, нежели простое отношение величин
В таком представлении итог с отклонением в десятом знаке будет выглядеть примерно как 0,9999999991, что ни при каком представлении не позволяет уменьшить разрядность записи числа.
Обнаруженное совпадение значений означает, что полное значение интеграла можно получить, вычислив его на любом интервале в 25% и умножив результат на 4. Для большей наглядности строим пять графиков: k(n), s25...s100. Поскольку значения колонок таблицы практически тождественны, все графики слились в одну линию.
Рис.5.7. Независимо от дискретности численного интеграла (5.8) все его интервальные значения равны друг другу
В связи с полученным результатом возникает вопрос: а если взять интервалы в 10% от общего? Проверка показала, что и в этом случае значения в каждой из соответствующих колонок даёт кратную, десятую часть от общей величины интеграла. Кроме того, при разных значениях нижнего предела и неизменном значении интервала Δ, то есть, превышения верхнего предела над нижним, значение интеграла в уравнении (5.7) приблизительно равно одной и той же величине
Эти значения интеграла лишь незначительно отличаются друг от друга на всех пределах и коротком интервале Δ. При исследовании на малых значениях дискретности замечено, что график интервальных интегралов имеет слегка выпуклую вверх форму, подобную рис.5.4. Указанный график - это зависимость значения интеграла для некоторого интервала Δ от величины нижнего предела. Высота этой выпуклости составляет практически доли процента от максимального значения.
Естественно, что такой же вопрос возникает и в отношении предельно короткого интервала - 1%. Интервалы целочисленные, поэтому мы не можем взять участок, интервал, меньший единицы. И вновь видим: значение интеграла на этом участке кратно равно значению интеграла на всём интервале. Если задать Δ = 0, то есть, единичный интервал, и точку интегрирования n/4, то уравнение интеграла с подстановочными и групповыми параметрами довольно точно описывается выражением общего вида
Обращаем внимание: оба предела внешнего интеграла равны n/4, то есть, интервал интегрирования равен единице. Здесь уместно вспомнить, что ранее мы такую возможность уже рассматривали и пришли к неутешительному выводу: полученный итог был явно неверным. Почему же сейчас мы вновь рассматриваем этот вариант и находим его полезным? Дело в том, что ранее мы нашли решение в предельном виде, сведя к единичным интервалам все интегралы. Здесь же мы опытным путём обнаружили, что третий, внешний интеграл допускает такое упрощение сам по себе, не ведя к противоречию, по крайней мере, на этом этапе.
Ещё раз напомним, что на этом этапе вычисления интеграла главной задачей является определение такого значения дискретности, превышение которой не приводит к его существенному уменьшению. Следует согласиться с тем, что при указанных в (5.9) условиях, пределах интегрирования значения интегралов kn являются допустимыми для сравнения. После определения предельного значения дискретности можно вычислить точное значение интеграла (5.4), хотя, видимо, и значение (5.9) можно рассматривать как приемлемо точное. Результаты вычислений (5.9) для разных дискретностей приведены в следующей таблице таб.5.7.
Таб.5.7. Вычисленные значения числового интеграла (5.9) уменьшаются, а его тенденция асимптотически стремится к единице.
Первая колонка - это значение дискретности, уменьшенное для компактности записи в 100 раз. Значения выбраны в геометрической прогрессии: каждое следующее в 2 раза больше предыдущего. Вторая колонка - вычисленный интеграл kn, а третья - отношение предыдущего значения интеграла к текущему, то есть, величина, которую мы назвали тенденцией интеграла. Для наглядности, данные из таблицы таб.5.7 приведены в графическом виде, в виде диаграммы рис.5.8.
Рис.5.8. Графики интеграла и тенденции, полученные числовым интегрированием выражения (5.9). Графики асимптотически стремятся к минимальным значениям.
На диаграмме видно, что график тенденции изменяется довольно медленно, но явно стремится к асимптоте - единице. Вычисление интеграла в двух последних строка таблицы таб.5.7 длилось, соответственно, почти 9 и 35 часов. Ожидается, что для вычисления следующего интеграла с дискретностью 8'192 потребует неприемлемо большое время - 140 часов.
Казалось бы, ситуация близка к тупиковой. Достигнутая погрешность превышает 7%, что мы считаем недостаточной точностью. Даже если и удастся дождаться вычисления 140-часового, недельного интеграла, то уже следующий потребует почти месяца вычислений. Вряд ли задача заслуживает таких больших затрат времени. Однако вновь проблема нашла неожиданное, случайное решение. Выше мы обнаружили, что полная величина интеграла может быть сформирована на отдельном участке, например, в 25% от общего предела интегрирования - таблица таб.5.5.
Предположительно это ведёт к сокращению времени интегрирования в 4 раза. Последние 4 строки таблицы таб.5.7 были повторно вычислены именно этим способом. Полученные ранее детальные вычисления показали, что все единичные интервалы практически тождественны - рис.5.3, таб.5.2.
Для контроля и наглядности вычисление интегралов в нашем программном комплексе мы сопровождаем выводом на экран текущего значения переменной и времени, прошедшего от начала интегрирования. Для ещё большей наглядности, а по сути, просто для красоты, выводим и текущее состояние диаграммы. Например, в процессе заполнения таблицы таб.5.7 на дискретности 512 (полное, точное значение 51'200. Напомним, что здесь и в некоторых других местах для компактности записей мы используем укороченную запись дискретностей, без двух нулей в конце) почти в самом начале вычисления интеграла (5.6) диаграммы выглядела, как показано на рисунке
Рис.5.9. Начальный, растущий этап численного вычисления интеграла (5.9), до появления скачка
Интеграл монотонно возрастает, а тенденция столь же монотонно уменьшается, стремясь к единице. При этом значения шкал, осей синхронно, автоматически меняются в соответствии со значениями крайних правых точек графиков. Наблюдение велось эпизодически, от случая к случаю, вплоть до завершения вычислений. Никаких особых изменений не ожидалось. Но произошло традиционное "вдруг". Событие выглядело весьма нелогично, поэтому в журнале исследований появилась запись:
"Совершенно случайно пойман скачок значений тенденции: в течение 6 минут шкала изменилась с 400 до 4. Это не очень хорошее обстоятельство. Видимо, эти вычисления следует перепроверить".
Действительно, при очередном, назовём его первым, наблюдении график выглядел, как показано на рис.5.9, но через 6 минут он кардинально изменился и стал выглядеть следующим образом
Рис.5.10. В момент прохождения точки скачка интеграл (5.10) и тенденция скачкообразно меняют значения, которые в дальнейшем практически не меняются.
Сомнения в его корректности понятны: графики рис.5.9 и рис.5.10 не имеют практически ничего общего, кроме названий. Однако оба они получены при вычислении уравнений, сформированных из уравнения (5.9):
По сути, внешний интеграл в уравнении теперь присутствует номинально, как разовое вычисление. Немного изменим алгоритм вычислений, для чего перепишем уравнение (5.9) в следующем виде:
Смысл этого нового интеграла в том, что мы вычисляем его при некотором разовом значении Δ внешней переменной, но верхний предел следующего интеграла равен этому значению плюс некоторое дополнительное число шагов, на котором образуется установившееся значение интеграла. Другими словами, мы вычисляем интеграл (5.10) с пределами внешнего интеграла Δ и прерываем вычисление, когда переменная во втором интеграле достигнет значения Δ + δ. Продолжение вычислений уже не требуется, поскольку вследствие обнаруженного скачка рост значения интеграла практически прекращается.
Приведённые рассуждения в значительной степени являются предположениями. Поэтому следует их проверить повторными вычислениями. Естественно, проверка на этом интервале дала тот же результат, графики. На других интервалах с этой же компактной дискретностью n = 512 получены схожие графики, значения. Вычисления на разных интервалах и дискретностях с подобными условиями сведены в таблицу
Таб.5.8. Таблица значений интеграла kn для разных дискретностей n = 400...204'800 и заданных точек скачка Δ = n%.
Левая, самая первая колонка таблиц n% - это значения точек скачка, значения Δ. Шапки таблиц n4...n2048 - это дискретности, под которыми в колонке приведены значения интегралов (5.10), и соответствующие им погрешности D4...D2048. Погрешности определяем, используя приведённый выше алгоритм, как отношение значений интегралов соответствующих скачков к максимальному значению интеграла для скачка 91%:
Поскольку в интеграле (5.10) мы использовали "плавающее" значение параметра δ, равное 10% от величины скачка, последнее значение скачка Δ мы выбрали равным 91%, иначе интервал интегрирования для последнего значения скачка окажется короче остальных. В таблице таб.5.8 видим, что независимо от этих пределов, от величины скачка, значения интеграла и тенденций практически одни и те же. Эффект скачка обнаружен на дискретности 512, поэтому для обобщения приведём диаграмму рис.5.11 в виде компиляции фрагментов графиков вблизи значения пределов интегрирования, скачка.
Рис.5.11. Компиляция диаграмм значений численного интегрирования выражения kn для дискретности n = 51'200 и заданных точек скачка Δ.
На рисунке рис.5.11 числа с процентом вверху - это значения пределов на горизонтальной шкале, то есть, шкала - это текущие значения переменной второго, среднего интеграла (5.6) в процентах от дискретности. Интеграл и тенденция плавно меняются вплоть до точки верхних пределов второго интеграла, после чего скачком меняют свои значения и дальше практически остаются неизменными. Шкала значений интеграла, правая вертикальная оказывается настолько сжатой, что его плавный рост до точки скачка практически выродился в горизонтальную линию.
Первый кадр, 2% - это вычисление по второй переменной - переменной y - от начала интервала до некоторого значения после точки скачка, когда значения графиков стабилизировались, стали меняться незначительно. Последний кадр, 75%, соответственно, заключительный участок от точки скачка до завершения вычислений. Очевидно, что такую заключительную "полку" имеют графики и на остальных кадрах.
В той же таблице таб.5.8 в дополнение к дискретности 512 проведены результаты ещё нескольких тестов, теперь уже с другими дискретностями. Например, график для дискретности 2'048 (полное значение - 204'800) выглядит следующим образом
Рис.5.12. Процесс численного интегрирования выражения k2'048 для дискретности 204'800 и c заданным значением скачка Δ = 2'048. В момент прохождения точки скачка интеграл и тенденция скачкообразно меняют значения, которые в дальнейшем практически не меняются.
Точка скачка Δ на диаграмме выбрана равной 1%. Все признаки скачка налицо: после прохождения его точки графики стабилизировались на новых значениях. Объяснить возникновение таких скачков, признаемся, оказалось задачей трудно разрешимой. Нерешённой. Видимо, это специфическое свойство сформированных интегралов. Но польза от него очевидная. Использование точки скачка позволяет вычислить интеграл за существенно меньшее время, чем непрерывное интегрирование. Видим, что интеграл растёт после точки скачка буквально всего лишь на протяжении нескольких шагов. Тем не менее, даже короткого шага после скачка достаточно для сравнения последовательных дискретностей. Точность собственно интегралов и тенденций в этом случае не играет решающей роли: необходимо найти самую большую дискретность приемлемого схождения интегралов. Для реализации этого режима интеграл и был преобразован к выражению (5.10).
Для обоснования возможности использования интеграла (5.10), оценим погрешности, возникшие при переходе от исходного четырёхкратного интеграла (5.4) к преобразованному двойному интегралу (5.10). С этой целью ещё раз обратимся к таблице таб.5.8. Сначала отметим, что по тем или иным соображениям в процессе вычислений мы выводили значения интегралов в разных масштабах. В таблице сравнения таб.5.8 мы учли эти масштабы. В качестве эталонных значений, с которыми сравниваем все остальные, мы используем самое первое значение, полученное при вычислении наименее преобразованного, наименее упрощённого интеграла. Например, для дискретностей 100 - 800 эти значения представлены в таблице таб.5.1. Все данные в таблице приведены в масштабе 106, использованном в большинстве вычислений.
Поскольку далее мы планируем вычислять интегралы со скачком в 1%, рассмотрим эту строку в таблице таб.5.8. Очевидно, со скачком в 91% значение интеграла и тенденции наиболее точное, судя по погрешностям. Чем меньше величина скачка n%, тем сильнее это значение отличается от значения скачка 91%, которое можно назвать базовым, поскольку оно просто наиболее близко к полному интервалу интеграла. Для малых дискретностей базовым значением являются значения, вычисленные на полных интервалах, без скачка. Видим в таблице, что погрешности в каждой колонке дискретностей растут снизу вверх. Но нас интересует рост погрешностей при скачке в 1% слева направо. Видим, что чем выше дискретность, тем выше погрешность, отклонение значения интеграла от его базовой величины. Конечно, желательно иметь колонки скачков и для дискретностей значения большего, чем 2048. Однако это неприемлемо длительные вычисления. Поэтому рассмотрим своеобразную тенденцию погрешностей для скачка в 1% в зависимости от представленных в таблице дискретностей. Диаграмма этой зависимости показана на следующем рисунке рис.5.13.
Рис.5.13. При увеличении дискретности рост погрешности замедляется
Видим, что погрешность растёт, но растёт с замедлением. Это значит, что на самых больших дискретностях погрешность можно ожидать в приемлемых границах, вероятно, не более 1%. В таблице таб.5.8 эта погрешность не превышает 0,3%. Отметим также, что сравниваем мы наихудшие погрешности, отношение интеграла на минимальном, самом грубом скачке, равном 1%, к интегралу, дающему наибольшую точность на предельно большом скачке в 91%. На следующем рисунке рис.5.14 это хорошо видно: с ростом величины скачка погрешность падает, но и наихудшие из них на скачке 1% имеют достаточно малые значения.
Рис.5.14. При любой дискретности погрешность растёт с уменьшением величины скачка
Из этого можно сделать вывод о корректности, допустимости вычисления интегралов именно на этом минимально возможном скачке. Если выявленное в итоге расхождение интегралов даже с грубыми значениями будет малым, то полное значение интеграла, очевидно, тем более будет входить в границы этой погрешности. На основе рассмотренного подхода, обоснований получен набор значений интегралов и тенденций для разных дискретностей с использованием эффекта скачка, равного 1%. Данные сведены в таблицу таб.5.9. За приемлемое время получены интегралы с дискретностями вплоть до 3'276'800. Заметим, что значение параметров с этой максимальной дискретностью получено за время почти 13 часов.
Таб.5.9. Таблица значений коэффициента k(d), интеграла (5.10) в зависимости от скачка
Как видим в таблице, тенденция на дискретности 3'276'800 уменьшилась до приемлемого значения: два смежных интеграла различаются менее чем на 1%. Следовательно, учитывая эти значения тенденций, за окончательное значение интеграла можно принять его значение с этой дискретностью, пусть и с некоторой погрешностью. На основе данных таблицы таб.5.9 на рис.5.15 приведены графики интеграла и тенденции как зависимости от дискретностей.
Рис.5.15. Итог численного интегрирования выражения (5.10) для всего диапазона дискретностей и со значением скачка Δ = n/100. Интеграл достиг своего предельного, минимального значения.
Таблица таб.5.9 и диаграммы рис.5.15 явно свидетельствуют о корректности полученных данных. Следует признать, что последняя строка таблицы приемлемо соответствует ожидаемым точным, реальным значениям интеграла и тенденции. Вместе с тем, можно возразить, что эти значения предельными всё-таки не являются. Очевидно, что интеграл при увеличении дискретности будет и дальше уменьшаться. Но до какой величины? Оценить её значение можно с помощью следующего приёма. Введём такие понятия как граница функции Fn и функциональная асимптота fn. Границей функции мы назовём такую функцию, значения которой и поведение которой нам известны наперёд с максимальной точностью и, кроме этого, функция по всем значениям близка к исследуемой. Удобнее всего, если она всегда будет меньше неё, будет границей снизу. Соответственно, функциональной асимптотой fn мы назовём функцию, равную разнице значений исследуемой функции, здесь это функция k(d), и границы функции Fn. Очевидно также, что для нашей исследуемой функции, для интеграла (5.10) удобнее всего использовать в качестве границы функции показательную функцию с основанием 1/2, поскольку мы выбрали в качестве шкалы дискретностей как раз эту последовательность. На рис.5.16 приведены графики этих трёх функций.
Рис.5.16. Использование границы функции и функциональной асимптоты
Использование логарифмической с основанием 2 шкалы дискретностей, горизонтальной оси позволило показать на этой оси реальные значения исследованных дискретностей (напомним, с отброшенными двумя нулями в конце). Мы отрегулировали границу функции Fn таким образом, чтобы у неё была асимптота fn, лишь немного ниже видимой, предполагаемой асимптоты графика интеграла k(d). Показательная функция имеет чётко определённую асимптоту, предел, который мы декларативно установили равным 0,0002, то есть, примерно на 0,0001 меньше, чем предполагаемая асимптота графика интеграла. Использованные в данном примере функции имеют вид
Используя указанный алгоритм построения функциональной асимптоты, строим и её график. В пределах визуальной точности, в пределах точности используемого программного комплекса мы видим, что график функциональной асимптоты выродился в горизонтальную линию. В абсолютных значениях отклонение асимптоты от минимального значения на пяти последних интервалах уменьшилось с 3% до 0,003%. Это означает, что график интеграла асимптотически стремится к определённому значению, равному 0,0002 + 0,0001 = 0,0003. Следовательно, это значение и следует принять как окончательное значение интеграла. Напомним, что значение вычисленного интеграла - это коэффициент kh. Этот коэффициент, равный с учётом использованного масштаба 3·10-4 : 106 = 3·10-10, означает, что сила гравитационного притяжения двух точек с массами, равными массам пластин, в 30'000'000'000 раз превышает таковую между пластинами, зеркалами.
Поясним использование термина "зеркала". Рассмотренная задача инициирована изучением эффекта Казимира, описывающего притяжение зеркал под действием, как считается, виртуальных частиц. Уравнение Казимира описывает относительную силу, то есть, силу, зависящую от площади зеркал. Эта сила обратно пропорциональна четвёртой степени расстояния между зеркалами. Однако и относительная сила гравитационного притяжения Ньютона в этом же варианте также обратно пропорциональна четвёртой степени расстояния между плоскими поверхностями. То есть, при таком подходе уравнение Ньютона выглядит как уравнение Казимира, отличаясь от него только константой
Как видим, формально два уравнения выглядят тождественно, поэтому мы и использовали в задаче термин "зеркала". В уравнении закона Ньютона площадь S взята равной квадрату расстояния между зеркалами, а коэффициент kh - полученное значение интеграла. В наших интегральных расчётах мы использовали площадь зеркал - 100 мм2, а расстояние между ними - 0,001 мм. Вычисленный коэффициент kh мы используем по той причине, что уравнение Ньютона описывает взаимодействие точечных объектов. Притяжение зеркал, как мы определили, можно выразить через взаимодействие точек той же массы.
В приложении П1 в демонстрационных целях приведён код процедуры числового интегрирования уравнения (5.10).
6. Притяжение точки к шару
Известно и зачастую подчёркивается явно, что закон всемирного тяготения Ньютона справедлив только для точечных тел или тел, расстояние между которыми многократно превышает их размеры. Рассмотрим, чему равна сила притяжения к шару точки, находящейся на некотором удалении от поверхности шара, то есть, за пределами применимости закона Ньютона в классической формулировке.
Рис.6.1. Притяжение к шару точки, находящейся на удалении от его поверхности, равном радиусу шара
Если бы вся масса шара была сосредоточена в его центре, то сила гравитационного притяжения между двумя одинаковыми массами m в случае рис.6.1 определялась бы уравнением Ньютона
Однако масса шара равномерно распределена по его объёму, поэтому для определения силы притяжения внешней точки m к такому шару нужно найти сумму её притяжения ко всем элементарным точкам dm шара. Расстояние между точками в каждой паре определяется координатами точки dm внутри шара. Проведём через произвольную точку 0' шара плоскость, ортогональную к координатной оси x. Эта плоскость выделит в шаре круг некоторого радиуса rmax. На этом круге, используя полярные координаты, построим окружность, обруч радиуса r. Расстояние R между любой точкой dm обруча, имеющей полярный угол φ, и внешней точкой m определяется соотношением
В общем случае радиус-вектор r произвольного обруча на круге изменяется в диапазоне от 0 до rmax. Полярный угол φ меняется от нуля до 2π. Максимальное, предельное значение радиус-вектора равно rmax
Объём элементарной массы dm на обруче равен
Обозначив через q плотность вещества шара, величина которой определяется массой шара и его объёмом, находим массу точки dm
C учётом того, что радиус R образует с осью x некоторый угол α, горизонтальная составляющая силы гравитационного притяжения масс dm и m вдоль оси x равна
Угол α, его косинус определяется соотношением
Подставляем dm и косинус в выражение для дифференциала силы
Результирующая сила притяжения точки к шару равна сумме всех элементарных сил, то есть, интегралу от полученного выражения
Разделим интегралы по переменным
Внутренний интеграл по переменной φ можем взять сразу же, поскольку подынтегральная функция от этой переменной не зависит. Заметим, что по этой переменной φ мы не меняем нижний предел с нуля на единицу, поскольку в данном случае вычисляем его как интеграл математический, а не численный, компьютерный.
Добавим множитель, чтобы выделить гравитационную силу между двумя точками
Разделим интеграл на два сомножителя
Первый сомножитель - это сила гравитационного притяжения FH между двумя точками массой m (6.1)
Назовём безразмерный второй сомножитель коэффициентом kh
Для наглядности развернём интегралы
Вычислим этот интеграл методом численного интегрирования. Для этого рассмотрим ряд дискретностей, то есть, разбиения переменных на элементарные участки n = 1000, 2000, 3000... Приведём уравнение (6.2) к виду уже сформированных процедур Excel. Для этого заменим переменные. Во внешнем интеграле используем новую переменную y
Заметим, что и в данном случае мы не меняем нижний предел по переменной y с нуля на единицу, поскольку в данном случае при интегрировании должны использоваться оба крайних значения: y = 0, соответствующее значению x = -R0, и y = n, соответствующее значению x = +R0. Верхний предел внутреннего интеграла становится равным
С этой заменой переменных интеграл (6.2) принимает вид
Теперь меняем переменную r во внутреннем интеграле на освободившуюся переменную x
С учётом этой замены интеграл принимает новый вид
Сделаем последовательно несколько преобразований
Введём в рассмотрение параметр p = h/R0
В таком виде уравнение является недостаточно наглядным, поэтому перепишем его в более наглядном и компактном виде
Для удобства формирования компьютерной процедуры собираем вместе уравнение интеграла и все его подстановочные групповые параметры
Выведенное уравнение позволяет рассмотреть частные варианты, специфические соотношения между радиусом шара и расстоянием точки над его поверхностью. Очевидно, наиболее интересными вариантами являются три: точка лежит на поверхности шара; точка находится над поверхностью шара на высоте, равной радиусу шара; и точка находится на достаточно большом удалении от поверхности шара. Рассмотрим эти варианты. Сначала рассмотрим "равновесный" вариант h = R0. Данные численного интегрирования сводим в таблицу
Таб.6.1. Притяжение к шару точки, находящейся над его поверхностью на высоте, равной радиусу шара
Видим, что в таблице коэффициент сравнения сил - интеграл k(n) стремится к единице при росте дискретности n интеграла. Тенденция также стремится к единице, подтверждая, что с увеличением дискретности n отношение сил стремится к единице. Гравитация шара эквивалентна гравитации точки в его центре с той же массой. По данным таблицы строим диаграммы тенденции и интеграла.
Рис.6.2. Диаграмма тенденции и коэффициента силы притяжения к шару точки, находящейся над его поверхностью на высоте, равной его радиусу
Теперь рассмотрим весьма интересный вариант: точка находится непосредственно на поверхности шара. Нетрудно догадаться, что эта ситуация напоминает наблюдателя, находящегося на поверхности Земли. Для нулевого удаления точки от сферы необходимо внести в уравнение интеграла некоторые исправления, устраняющие деление на ноль. Но мы поступим иначе, просто сделаем высоту предельно малой: h = 0,00001R0.
Таб.6.2. Притяжение к шару точки, лежащей на его поверхности
Исходя из смысла понятия Интеграл в таблице как коэффициента сравнения сил, и данном случае видим, что сила гравитационного притяжения точки к шару равна силе её притяжения к точке в центре шара с его массой. Интеграл k(n) стремится к единице при росте дискретности интеграла. Тенденция также стремится к единице, подтверждая, что с увеличением дискретности n отношение сил стремится к единице. Другими словами, гравитация шара эквивалентна гравитации точки в его центре с той же массой. По данным этой таблицы строим новые диаграммы тенденции и интеграла
Рис.6.3. Диаграмма тенденции и силы притяжение к шару точки, лежащей на его поверхности
Диаграмма подтверждает сделанный выше вывод: гравитация шара эквивалентна гравитации точки в его центре с той же массой. Таблица подтверждает, что сила гравитационного притяжения точки к шару равна силе её притяжения к точке в центре шара с его массой. Коэффициент сравнения сил - интеграл k(n) стремится к единице при росте дискретности n интеграла.
Отметим, что в первом варианте при равенстве h = R0 графики стремятся к асимптотам сверху, а во втором, h = 0 - снизу. Последний вариант является обобщённым, для произвольно выбранного значения высоты точки над поверхностью шара. Точка удалена от поверхности сферы на большое расстояние, h = 5R0. Полученные данные вносим в новую таблицу и строим по ним диаграмму рис.6.4
Таб.6.3. Притяжение к шару точки, находящейся над его поверхностью на высоте, равной пяти радиусам шара
Рис.6.4. Диаграмма тенденции и коэффициента силы притяжение к шару точки, находящейся над его поверхностью на высоте, равной пяти его радиусам
Хотя мы повторили это уже несколько раз в отношении других диаграмм и таблиц, но и в этом случае вновь будет уместно повторить: гравитация шара эквивалентна гравитации точки в его центре с той же массой.
Используя числовое интегрирование, мы получили три частных решения, каждое из которых отчётливо показало: гравитационное притяжение точки к шару таково, будто точка притягивается к центру шара, к условной точке, в которой сосредоточена вся его масса. Это следует из того, что к единичной асимптоте стремится не только тенденция, но и коэффициент. Коэффициент - это отношение силы притяжения между шаром и точкой к силе притяжения между точкой и центром шара, в котором сосредоточена вся его масса.
Следует признать, что такой результат выглядит весьма странно, непривычно, если вспомнить ограничения, накладываемые на закон Всемирного тяготения относительно размеров взаимодействующих тел. Поэтому можно назвать удачей возможность проверить это решение аналитически. Действительно, интеграл (6.2) можно свести к табличным, то есть, получить решение в общем аналитическом виде. Поскольку ожидаемое аналитическое решение может вызвать недоверие, приведём его в предельно детальном виде. Внутренний интеграл по переменной r к табличному виду сводится довольно просто.
Подставляем пределы и меняем знак, меняя местами слагаемые
Подставляем найденное решение в исходное уравнение
Разделяем уравнение на два интеграла
Решаем их поочерёдно
Решаем второй интеграл
Этот интеграл также разбиваем на два
Также решаем их по очереди
Этот интеграл вновь решаем последовательными преобразованиями
Раскрываем скобки под корнями
Собираем подобные члены
Собираем оставшиеся подобные члены
Упрощаем выражение
Собираем коэффициент kh
Теперь вычисляем последний интеграл
Интеграл является табличным вида
Подставляем в табличную форму величины нашего интеграла пока без пределов интегрирования
Делаем простейшие преобразования, упрощения
Подставляем пределы интегрирования и преобразуем
Находим значение коэффициента kh4 и собираем окончательное значение коэффициента kh
Приведём первичное выражение для сравнения сил притяжения точки и шара
Полученное соотношение означает, что независимо от удалённости от шара, точка всегда притягивается к его центру. Закон Всемирного тяготения Ньютона в классической формулировке в точности описывает притяжение точки к однородному объекту сферической формы, шару
Решение получено в общем виде, поскольку для численного интегрирования мы не использовали никаких конкретных числовых значений радиуса шара и расстояния между точкой и шаром.
Притяжение двух шаров
Ещё раз вспомним известное мнение, что закон всемирного тяготения Ньютона справедлив только для точечных тел или тел, расстояние между которыми многократно превышает их размеры. Однако приведённый только что точный анализ показал, что притяжение точки к массивному шару эквивалентно её притяжению к такой же точке с массой шара, находящейся в его центре. Проведём дальнейший анализ и рассмотрим, чему равна сила притяжения двух шаров друг к другу. Воспользуемся тем же рисунком, что и при исследовании притяжения точки к шару, но рассматривать его будет условно и зеркально, считая, что рассмотренный выше шар R0 теперь находится справа и обозначен как шар радиуса R1.
Рис.6.5. Притяжение к шару справа точки, являющейся центром другого шара (слева) и находящейся на некотором удалении от его поверхности
В задаче точка-шар мы пришли к выводу, что любая точка притягивается к шару, который можно рассматривать как такую же точку. На рисунке рис.6.5 точку m справа, находящуюся в начале координат xz", будем рассматривать как массивный центр такого шара R1, подобного тому, что рассмотрен ранее, в предыдущем разделе. Считаем, что вся масса этого шара R1 сосредоточена в этой точке, в его центре. Очевидно, радиус R1 этого шара может быть любым, поскольку в уравнение притяжения к нему любой внешней, в том числе любой точки шара R0, который мы в данном случае считаем новым, этот радиус входит как часть расстояния между центром шара R1 и произвольной точкой, точкой шара R0. Независимо от удалённости от поверхности этого шара R1, внешняя точка, любая точка шара R0 всегда притягивается к его центру
Здесь, как показано на рисунке, R - это расстояние между элементарной точкой нового шара R0 и центром шара R1. Мы специально привели исходный рисунок в том же виде, как и в исходной задаче шара и точки. Обозначения нового шара R0 соответствуют обозначениям шара в предыдущем разделе для совместимости выкладок.
Очевидно, все приведённые в предыдущем разделе рассуждения следует повторить слово в слово, поскольку здесь мы опять имеем ту же саму ситуацию: притяжение точки - центра некоторого шара (теперь он справа) и произвольной точки, теперь это элементарные точки нового шара слева. Повторим, шар справа показан иллюстративно, его радиус может быть любым, меньшим h, поскольку притягиваемой гравитационной точкой m является его центр.
Очевидно, что все рассуждения для пары шар-точка приведут к такому же заключению и, как следствие, к выводу: сила притяжения двух однородных шаров зависит только от их масс и расстояния между их центрами. Размеры шаров, их радиусы могут быть любыми, вплоть до соприкосновения шаров. Это означает ошибочность представления о том, что закон Всемирного тяготения Ньютона справедлив только для точечных тел или тел, размеры которых значительно меньше расстояния между ними. Закон в классической формулировке справедлив для тел шарообразной формы при любом расстоянии между ними.
Следует отметить: этот один и тот же результат получен аналитически и численным интегрированием, то есть, двумя принципиально разными способами, которые явно показали, что притяжение двух шаров не зависит от их размера, а определяется исключительно их массами и расстоянием между их центрами. Эти два результата не просто повторяют друг друга, они взаимно утверждают истинность этих результатов, выводов.
7. Притяжение двух дисков
Исследуя гравитационное притяжение сферических объектов, мы обнаружили довольно странное обстоятельство. Вопреки существующему мнению, два шара притягиваются с такой же силой, как и две точки, находящиеся в центрах этих шаров и имеющих их массу. Возникает естественный вопрос: есть ли другие объекты с подобными свойствами? Например, с какой силой притягиваются два диска, отдалённые "родственники" шаров? Чтобы выяснить этот вопрос, рассмотрим схему, приведённую на следующем рисунке.
Рис.7.1. Расчёт силы гравитационного притяжения двух дисков
Два бесконечно тонких диска одинакового диаметра находятся на некотором расстоянии d друг от друга. Сила их гравитационного притяжения равна сумме притяжения всех элементарных точек одного диска ко всем элементарным точкам другого диска. Каждая такая элементарная, дифференциальная сила между двумя произвольными точками равна
Дифференциалы масс точек определяются уравнениями
Дробь, отношение массы диска к его площади - это поверхностная плотность диска. Точки на каждом диске рассматриваем в полярных координатах, что в данном случае более удобно. Расстояние между этими произвольными элементарными точками равно
Здесь мы использовали закон косинуса для треугольника, образованного радиус-вектором нижней точки и проекцией на нижний диск радиус-вектора верхней. Подставляем эти величины в уравнение (7.1)
Эта дифференциальная сила направлена вдоль радиуса, соединяющего точки, а нас интересует сила, направленная ортогонально к дискам. На рисунке видим, что отношение этих двух сил описывается пропорцией
Используя эту пропорцию, находим дифференциал силы, ортогональной к дискам
Суммарную силу гравитационного притяжения дисков находим интегрированием
Разнесём переменные
Анализируя рисунок, замечаем, что точка с углом α = 0 повторяется 2π раз, поэтому внутренний интеграл по углу α вычисляем сразу же
Упростим знаменатель в подынтегральной функции
Подставляем его обратно в интеграл
Как и раньше, нас интересует не собственно сила притяжения дисков друг к другу, а отличие этой силы от силы притяжения двух точек с массами дисков и на таком же удалении друг от друга. Для определения этого различия сил преобразуем интеграл, умножив его на единичный коэффициент - отношение квадрата расстояния d2 к самому себе
Выделяем в отдельный множитель величины, описывающие указанную силу FH притяжения двух точек друг к другу
Второй сомножитель с интегралами обозначим коэффициентом k
Если этот коэффициент при любых значениях параметров будет равен единице, то это означает, что два диска притягиваются так, будто являются точками. Для компьютерного вычисления интеграла и преобразования уравнения к формату уже сформированных процедур Excel сделаем ставшую традиционной замену переменных. Замену сделаем раздельно для внешнего интеграла по углу φ и для двух внутренних. Диапазоны изменения исходных и переменных для замены соответствуют друг другу. Замечаем, что по углу φ интеграл симметричный, поэтому угол берём в половину окружности - от нуля до 180 градусов, интеграл, соответственно, удваиваем
Связь между новыми и прежними переменными и их дифференциалами имеет вид
Уточним: в этих уравнениях m - это дискретность угла φ, а не масса. После подстановки уравнение (7.2) приобретает вид
В процессе пробных вычислений обнаружено появление чрезмерно больших величин, приводящих к сбою. Для их устранения преобразуем уравнение для процедуры. В теле, в коде процедуры оставляем дроби - отношение линейных переменных к дискретности, остальные дроби выносим за пределы интегралов
Сокращаем
Выделяем групповые параметры и подставляем их в уравнение
Итоги численного интегрирования (7.4) сводим в таблицу таб.7.1. Полученные данные показывают, что даже при низких, грубых дискретности угла m угла φ и дискретностях n радиусов x и y, значение интеграла одно и то же с достаточно высокой точностью.
Таб.7.1. Таблица численного интегрирования выражения (7.4)
Теперь рассмотрим, можно ли применять закон Всемирного тяготения Ньютона к конфигурации из двух таких тонких дисков. Используя минимальное, наихудшее значение дискретности угла m = 90 и дискретности радиусов обручей n = 400, вычислим интегралы для разных значений расстояния между дисками d = 1...20. Уравнение для численного интегрирования имеет вид (7.5). Добавлены подстановочные, групповые параметры.
Для сравнивания силы притяжения дисков с силой притяжения двух точек, добавляем уравнение тенденции, которое в данном случае показывает отношение силы притяжения дисков к силе притяжения двух точек той же массы.
Итог вычислений интеграла (7.5) приведён на рис.7.2. На диаграмме график k(d) показывает, что по мере роста расстояния d между дисками сила их притяжения приближается к силе притяжения Fн двух точек такой же массы: коэффициент k(d) асимптотически стремится к единице. Следовательно, закон Всемирного тяготения Ньютона для таких дисковых объектов справедлив только для расстояний между ними, существенно превышающих их поперечные размеры. На таких расстояниях объекты можно рассматривать как точечные. График FH на рис.7.2 демонстрирует обычное уменьшение силы гравитационного притяжения по мере роста расстояния между массивными объектами, в данном случае - точками.
Рис.7.2. Итоги численного вычисления интеграла (7.5).
Обратим внимание на интересную картину на рис.7.2. На начальном этапе, слева на максимально близком расстоянии между дисками, сила их притяжения прямо пропорциональна силе притяжения точек, хотя и меньше неё. Об этом свидетельствует также и коэффициент k(d) - на этом интервале он практически не меняется. Также весьма загадочно выглядит изгиб, провал на графике Fd - силе притяжения дисков. На некотором критическом расстоянии между дисками сила их притяжения достигает локального максимума и далее падает в том же темпе, что и сила притяжения точек. Это расстояние, после которого диски можно рассматривать как точки в классическом варианте закона Всемирного тяготения.
Отметим ещё одну особенность диаграммы. Сила притяжения по закону Ньютона обратно пропорциональна квадрату расстояний между взаимодействующими телами, но на рисунке рис.7.2 мы видим, что сила изменяется линейно. Почему? Такая линейная зависимость связана просто с выбранной нелинейной, логарифмической системой координат. Отбросив константы, силу притяжения можно условно описать уравнением
Вертикальная ось, шкала силы выбрана логарифмической. Следовательно, график силы имеет вид
На диаграмме ось интервалов также выбрана логарифмической
Подставляем в уравнение для y
Это график прямой линий, имеющей отрицательный наклон, то есть, направленной слева сверху вправо вниз, и проходящей через начало координат. Именно так и выглядит график силы Fн на диаграмме.
Ранее мы устанавливали нижние пределы интегралов в единицу, поскольку расчёты велись в отношении каждого дифференциального участка объекта - шара или отрезка. В данном случае вычисления мы ведём не по номеру какого-либо дифференциального объекта, а по его физическому положению. Радиусы в данной формулировке интеграла изменяются от нуля до R0 включительно. Переменные, заменившие исходные радиусы дисков мы обозначили параметром n. Следовательно, в данной версии интеграла мы обязаны сохранить пределы по радиусам полными - от нуля до максимального значения параметров - n. Напротив, угол α на половине диска изменяется строго от нуля до величины, меньшей 180 градусов на одну дискрету, то есть, мы обязаны вычислять интеграл, исключив эту лишнюю дискрету. Поскольку мы использовали подстановку переменных, то угол половины окружности стал измеряться в целочисленных дискретах, "штуках", которые мы обозначили переменной m. Следовательно, нижний предел мы должны увеличить на единицу. Либо на единицу уменьшить верхний. Мы использовали в (7.5) диапазон 0...(m - 1) для пределов внешнего интеграла. И, действительно, только в этом случае мы получаем корректное значение силы притяжения дисков: при их предельном сближении сила притяжения стремится к бесконечности, как и сила притяжения двух точек. Исключение из рассмотрения нулевых точек дискретностей n, приводит к сомнительному результату: при уменьшении расстояния между дисками сила их притяжения как бы уменьшается, причём даже как бы стремится к нулю. Понятно, что такой результат не может быть верным.
Тем не менее, отметим, что обнаруженная волна, изгиб графика Fd: падение силы притяжения, затем её кратковременный рост с последующим непрерывным падением, выглядит довольно странно. На диаграмме рис.7.2 с логарифмическими шкалами, крайний правый участок графика имеет наклон, стремящийся к 45 градусам. Это означает, что никакой асимптоты нет, и абсолютное значение силы стремится к нулю, становясь в пределе равной силе притяжения двух точек той же массы, так же стремящейся к нулю.
Точка в плоскости диска
Как видим, описание гравитационного взаимодействия двух дисков не совпало с законом взаимодействия двух точек, с законом Всемирного тяготения Ньютона для любых расстояний между ними. Можно предположить, что это относится исключительно к данной конструкции - параллельному расположению двух дисков. Поэтому рассмотрим другую конструкцию: диск и точка, подобную исходной схеме шара и точки.
Рис.7.3. Определение силы гравитационного притяжения точки к диску
Точка находится в плоскости диска на некотором удалении от него, за его пределами. Используем те же физические параметры и принципы. Сила гравитационного притяжения точки к диску равна сумме притяжения её ко всем элементарным точкам диска. Каждая такая элементарная, дифференциальная сила притяжения к произвольной точке dm диска равна
Дифференциалы масс точек диска определяются уравнениями
Дробь, отношение массы диска к его площади, как и выше - это поверхностная плотность диска. Точки на диске рассматриваем в полярных координатах. Расстояние между внешней точкой и произвольной точкой dm диска равно
Подставляем эти величины в уравнение (7.6)
Эта дифференциальная сила направлена вдоль радиуса R, соединяющего точки, а нас интересует сила, направленная вдоль оси x. На рисунке видим, что отношение этих двух сил описывается пропорцией
Используя эту пропорцию, находим дифференциал силы, параллельной оси x
Суммарную силу гравитационного притяжения внешней точки к диску находим интегрированием
Разнесём переменные по интегралам
Как и раньше, нас интересует не собственно сила притяжения диска к внешней точке, а отличие этой силы от силы притяжения двух точек с массами диска и внешней точки и на таком же удалении друг от друга. Для определения этого различия сил преобразуем интеграл, умножив его на единичный коэффициент
Преобразуем, выделив в отдельный множитель величины, описывающие указанную силу FH притяжения двух точек друг к другу
Второй сомножитель с интегралами обозначим коэффициентом kx
Если этот коэффициент при любых значениях параметров будет равен единице, то это будет означать, что точка притягивается к диску так, будто тот являются точкой с его массой, расположенной в его центре. Для компьютерного вычисления интеграла и преобразования уравнения к формату уже сформированных процедур Excel сделаем ставшую традиционной замену переменных. Замену сделаем раздельно для внешнего интеграла по углу α и для внутреннего интеграла. Диапазоны изменения исходных и переменных для замены соответствуют друг другу. Замечаем, что по углу α интеграл симметричный, поэтому угол берём в половину окружности - от нуля до 180 градусов, интеграл, соответственно, удваиваем. Связь между новыми и прежними переменными и их дифференциалами имеет вид
Вновь уточним: в этих уравнениях m - это дискретность угла α, а не масса диска или точки. После подстановки и замены нижних пределов, уравнение (7.7) приобретает вид
Далее мы делаем последовательно несколько преобразований. Преобразования достаточно простые, поэтому мы не приводим их здесь полностью.
В процессе преобразований замечаем повторяющийся фрагмент, который обозначаем переменной p
Заменяем в уравнении повторяющиеся фрагменты его обозначением, назначенной переменной
Добавляем к блоку уравнение для силы притяжения опорных точек и окончательно записываем
Вычисления будем производить для разных значений расстояний h > 0 между диском и точкой, и разных дискретностей n. Без нарушения общности примем R0 = 1. Вычисление интеграла (7.8) с выбранными параметрами и условно принятым Gm2 = 2 в уравнении для силы Fн показало практически один и тот же результат для разных значений дискретности
Рис.7.4. Сравнение численным интегрированием силы притяжения диск-точка и двух точек
Видим, что при увеличении удалённости h внешней точки от диска коэффициент kx = k(h) асимптотически стремится к единице, графику k(1). Это очевидно, поскольку на больших удалённостях диск можно рассматривать как точку, то есть, сила притяжения в точности равна силе притяжения двух точек, которую описывает график FH. При увеличении расстояния между диском и точкой эта сила падает вплоть до нуля. Таким образом, к этой системе закон Всемирного тяготения Ньютона применим лишь при значительном расстоянии между точкой и диском.
Точка над плоскостью диска
Мы пришли к выводу, что гравитационное взаимодействие диска и точки, лежащей в его плоскости, не описывается законом взаимодействия двух точек, законом Всемирного тяготения Ньютона. Однако есть ещё один вариант расположения точки относительно диска - над его плоскостью на осевой линии. Эта схема изображена на следующем рисунке
Рис.7.5. Притяжение к диску точки, находящейся над плоскостью диска на его оси
Если бы вся масса диска была сосредоточена в его центре, то сила гравитационного притяжения между двумя одинаковыми массами m в случае рис.7.5 определялась бы уравнением Ньютона
Однако масса диска равномерно распределена по его плоскости, поэтому для определения силы притяжения внешней точки m к такому диску нужно найти сумму её притяжения ко всем элементарным точкам dm диска. Расстояние между точками в каждой паре определяется координатами точки dm на поверхности диска. Рассмотрим на диске в полярных координатах круг, обруч некоторого радиуса r. Расстояние R между любой точкой dm обруча, имеющей полярный угол φ и внешней точкой m определяется соотношением
В общем случае радиус-вектор r произвольного обруча на круге изменяется в диапазоне от 0 до R0. Полярный угол φ меняется от нуля до 2π. Объём, вернее, элементарная площадь элементарной массы dm на обруче равна
Обозначив через q плотность вещества диска, величина которой определяется массой и площадью диска, находим массу точки dm
C учётом того, что радиус R образует с осью x некоторый угол α, горизонтальная составляющая силы гравитационного притяжения масс dm и m вдоль оси x равна
Угол α, его косинус определяется соотношением
Подставляем dm и косинус в выражение для дифференциала силы
Результирующая сила притяжения точки к диску равна сумме всех элементарных сил, то есть, интегралу от полученного выражения
Разделим интегралы по переменным
Внутренний интеграл по переменной φ можем взять сразу же, поскольку подынтегральная функция от этой переменной не зависит. Заметим, что и в данном случае по этой переменной φ мы не меняем нижний предел с нуля на единицу, поскольку вычисляем его как интеграл математический, а не численный, компьютерный.
Добавим множитель, чтобы выделить гравитационную силу между двумя точками
Сразу же выделяем первый сомножитель - это сила (7.9) гравитационного притяжения FH между двумя точками массой m
Второй, безразмерный сомножитель назовём коэффициентом kh
Прежде, чем строить график этой интегральной функции, заметим, что этот интеграл - табличный. Поэтому для сравнения дополнительно произведём его вычисление аналитически. Это позволит нам построить диаграммы, графики обоих решений. Решение интеграла начнём с замены интегральной переменной x = h2+r2
Возвращаем исходную переменную и добавляем пределы интегрирования
Данное аналитическое решения показывает, что отношение сил зависит и от размера диска и от расстояния между точкой и плоскостью диска. То есть, применение классической формы закона Ньютона, закона Всемирного тяготения в данном случае также невозможно. Использование этого закона возможно только при условии h >> R0. Это условие можно преобразовать в эквивалентное
Поскольку оба параметра, по сути, являются переменными, другой, эквивалентной записью этого выражения является его запись в виде предела, в котором отношение мы обозначим переменной p
Подставляем эту переменную в аналитическое решение
Закон Всемирного тяготения можно использовать только при выполнении условия (7.11). Соотношение (7.12) можно записать в завершённом виде
Другими словами, при росте отношения (7.11) значение коэффициента kh должно стремиться к единице
Однако напрямую доказать этот предел оказалось довольно сложно: в уравнении разница двух бесконечностей равна конечной величине, единице, что, по сути, означает неопределённость. Непосредственное вычисление предела даёт его нулевое значение, которое нас устроить не может.
Вместе с тем, два слагаемых в предыдущем уравнении определённо не равны друг другу, то есть, при любой величине p следует ожидать какое-то конечное значение разницы между ними.
В другом варианте, при прямом решении уравнения kh = 1 был получен довольно странный результат. Простые преобразования привели к вполне убедительному решению, конкретному конечному значению величины p, однако подстановка этого решения в уравнение не дало единичного результата (приложение П2). То есть, найденное решение не было решением уравнения. Функция оказалась с загадкой. Точное решение удалось получить, сделав довольно очевидные и простые преобразования. Перепишем уравнение (7.13)
Это равенство содержит в скобках две бесконечно большие величины, обозначенные как A? и B?, которые равны, тождественны, поскольку разница между ними равна нулю.
Поскольку слагаемые равны, предел их отношения при росте переменной должен стремиться к единице. Проверяем, так ли это, и убеждаемся, что это действительно так
Это решение является строгим, хотя поверхностный логический анализ может предложить иной исход. Действительно, уравнение (7.13) мы можем записать иначе:
Тогда преобразования приведут к пределу, имеющему несколько иной вид
Из чего, казалось бы, мы обязаны сделать вывод, что пределом kh при выполнении условия (7.11), стремлении p к бесконечности может быть любое число, то есть, предел по-прежнему неопределённый. Однако этот подход содержит очень тонкую, практически незаметную логическую ошибку. При вычислении как старого, так и нового предела мы отбрасываем конечную величину в сумме её с бесконечно большой величиной. Как правило, это допустимо, это называется отбрасыванием бесконечно малой величины. То есть, мы по умолчанию отождествляем это отбрасывание в двух формах записи пределов
Но в данном случае есть очень важная особенность. Фактически, отбрасывая бесконечно малые, во втором варианте предела мы отождествляем две величины - суммы в числителе и знаменателе, формально считая их равными
При поверхностном подходе это верно. Но какой бы большой ни была величина p2, это равенство при строгом логическом подходе является ошибочным. Действительно, перепишем его иначе
Из этого получаем
Но это явно неверное решение, поскольку коэффициент kh - число положительное. На самом деле исходное выражение для предела следует записать иначе, следующим образом
Теперь мы имеем полное право заменить тождественные величины в скобках на их эквивалентные тождественные значения
Именно это мы и сделали в исходном виде предела. Следовательно, полученное изначально решение является точным и строго логичным
Если в числителе вместо единицы записать произвольное число - a, то его отбрасывание как бесконечно малой становится логически недопустимым. Ещё одним весомым доказательством справедливости уравнения (7.13) можно назвать его прямое вычисление, вычисление kh. В пределах разрядности нашей системы, максимально возможных в ней чисел, решение (7.13) выполняется. В отношении рассмотренной системы это означает справедливость правила "закон Всемирного тяготения в классической формулировке выполняется, если размеры объектов много меньше расстояния между ними". Сила притяжения между диском и точкой над ним становится равной силе притяжения двух точек только при достаточно большом расстоянии между точкой и диском или когда размеры диска много меньше этого расстояния.
Теперь можно построить графики обоих решений. Для численного интегрирования воспользуемся косвенной заменой переменных. Будем измерять расстояния и радиус диска в предельно "мелких" единицах. То есть, назначим радиусу R0 любое достаточно большое значение, например, в ангстремах, что эквивалентно его замене на дискретность n большой величины:
Интегральная переменная r теперь определяется дискретностью радиуса R0. Величину h будем рассматривать как параметр интегрирования и диаграмму построим как функцию коэффициента kh(h). Поскольку эта зависимость для нас является основной, просто задаём некоторое достаточно большое значение радиуса и дискретности n = R0.
Рис.7.6. Сравнение относительной силы притяжения точки к диску
Как видим на диаграмме, с ростом удалённости точки от центра диска, относительная сила её притяжения к диску растёт до асимптотического предела, равного единице. В центре диска сила притяжения равна нулю, а по мере удаления от центра она приближается к силе притяжения точек той же массы. Закон Всемирного тяготения в классическом виде (7.9) для точки и диска соблюдается только при выполнении условия (7.11) - размеры диска много меньше расстояния между ним и точкой.
8. Притяжение двух стержней
Соосные стержни
Мы выяснили, что однородные массивные шары гравитационно притягиваются своими центрами, в которых как бы сосредоточена вся их масса. Это значит, что к таким объектам закон Всемирного тяготения Ньютона применим в его классической формулировке. Однако дискообразные тела под классическую формулировку закона уже не подпадают. Следует ожидать, это ограничение распространяется и на тонкие массивные стержни, осевые линии которых совпадают. Предполагаем, что сила гравитационного притяжения зависит от размеров стержней, но не определяется притяжением их центров. Рассмотрим схему гравитационного притяжения двух тонких стержней a и b, представленную на следующем рисунке.
Рис.8.1. Гравитационное притяжение двух стержней
Сила притяжения формируется из дифференциальных сил притяжения элементарных, дифференциальных отрезков стержней. Считаем, что диаметры стержней равны нулю, а масса дифференциальных участков определяется их длиной и линейной плотностью стержней. Сравним силу притяжения двух точек A и B стержней a и b, центры которых находятся на расстоянии h друг от друга. Массы стержней и эквивалентных им точек принимаем равными m. Базовая, сила притяжения двух точек по закону Ньютона равна
Силу притяжения стержней определяем как сумму дифференциальных сил притяжения между всеми точками одного стержня со всеми точками другого. Каждая такая элементарная, дифференциальная сила равна
Дифференциалы масс отрезков определяются дифференциалом их длины и линейной плотностью стержней, равной отношению массы стержня к его длине
Подставляем в исходное уравнение для dF
Расстояние между дифференциальными массами зависит от их координат, которые изменяются в следующих пределах
Исходя из этого, находим расстояние между любыми двумя дифференциальными массами
Подставляем в уравнение для дифференциала силы
Общую силу гравитационного притяжения стержней находим интегрированием
Разделяем интеграл на два по переменным
Для численного интегрирования присвоим константам определённые значения. Без нарушения общности установим Ra=Rb=R, при этом 0 < 2R < h. Равенство длин стержней расстоянию между их центрами пока исключаем, поскольку оно означает соприкосновение стержней и возможное появление в вычислениях бесконечностей. Также выбираем такие значения масс, чтобы обратить в единицу множитель
Подставляем этот множитель в уравнение
Теперь сделает традиционную замену переменных. Диапазоны изменения исходных и переменных для замены соответствуют друг другу.
Связь между новыми и прежними переменными и их дифференциалами имеет вид
Подставляем в уравнение новые переменные и сразу отбрасываем у них штрихи
Делаем преобразования, упрощающие вид интеграла
Окончательно записываем
Вновь обращаем внимание на диапазон пределов интегрирования. Мы выбрали абстрактную целочисленную дискретность n, которая означает, что диапазон определения функции разбит на n интервалов, участков. Но вновь напомним, что мы используем в VBA-процедурах конструкцию For...Next. Как мы отметили интуитивно и убедились в процессе вычислений, число шагов в цикле должно быть на единицу меньше. Действительно, при вычислении функции (8.3) мы получили завышенное значение интеграла. Простым сужением диапазона либо до 0...(n - 1), либо до 1...n решение приходило к естественному значению, поскольку просчитывались отрезки, существующие действительно. Очевидно, пределы 1...n в уравнении интеграла (8.3) выглядят более компактно, да и красиво.
Обращаем внимание: интеграл (8.2), судя по всему, является табличным, поэтому найдём его аналитическое решение для последующего сравнения с результатом численного интегрирования. Внутренний интеграл решаем последовательными преобразованиями
Подставляем в исходное уравнение и подставляем пределы
Преобразуем
Окончательно получаем интеграл, который также является табличным интегралом логарифмического вида
Подставляем наши параметры
Подставляем пределы и получаем итоговое уравнение
Окончательно записываем
Немного упростим запись этого решения
Для сравнения собираем в один блок три уравнения, добавив в него уравнения (8.1) и (8.3), исправив в последнем нижние пределы, как отмечено выше, на единичные
Строим пробные диаграммы этих трёх уравнений, для чего задаём значения параметров: h = 10. Переменный параметр R рассматриваем в максимальном диапазоне 0 < 2R < h.
Рис.8.2. Сравнения сил притяжения стержней и точек. Аналитический и интегральный графики
Видим, что даже выбранная непривычно низкая дискретность n = 200 даёт интегральное решение (красный график), идеально совпадающее с аналитическим (синий график). Зелёный график - это сила притяжения двух точек, график, которой фактически играет роль асимптоты, поскольку эта сила от длины стержней не зависит. Тем не менее, отметим, что на верхней границе интервала аналитический график немного расходится с численным. Связано это всё-таки с недостаточной дискретностью численного интегрирования. Предельное значение длины 2R = h ведёт к бесконечно большой силе притяжения.
Растущая сила притяжения двух соосных стержней означает, что их притяжение таково, будто они притягиваются не своими центрами тяжести, а некоторыми точками, которые в этом случае можно назвать гравитационными центрами стержней. Найдём положение этих центров, которые, очевидно, находятся к центру системы ближе, чем центры тяжести стержней. Рассмотрим два нижних уравнения из (8.4). Считаем, что сила притяжения стержней равна силе притяжения точек, которые находятся на некотором гравитационном удалении hg друг от друга. Это равенство запишем в следующем виде
Преобразуем равенство
Находим значение гравитационного радиуса
Окончательно запишем
Используя это уравнение, найдём аналитически некоторые предельные значения этого гравитационного радиуса.
Соприкосновение стержней: R = h/2
Если R = h/2, то есть, стержни касаются друг друга, гравитационный радиус стремится к нулю:
Это означает, что два соприкасающиеся стержня притягиваются так же, как и две точки той же массы, но находящиеся на нулевом удалении. Физически это ведёт к бесконечно большой силе притяжения.
Длина стержней стремится к нулю: R → 0
При стремлении радиуса R, длины стержней к нулю, аналитическое уравнение для гравитационного радиуса ведёт к неопределённости
Найдём предел этого отношения
Удобно найти предел подкоренного выражения
Предел отношения двух функций равен пределу отношения производных этих функций
Следовательно
Строго говоря, такой результат очевиден. При уменьшении длины стержней до нуля они превращаются в точки. Расстояние между точками по определению равно h, расстоянию между центрами стержней.
Произвольная длина стержней: R =0 ... h/2
Сначала рассмотрим два стержня, длины которых меньше расстояния между их центрами в 4 раза: h = 10, R = h/4 = 2,5
Два стержня длиной, равной четверти расстояния между их центрами, притягиваются как две точки, находящиеся на удалении друг от друга, равном примерно 9,3 их длин. Расстояние между центрами стержней равно 10, то есть, гравитационные центры стержней находятся друг к другу ближе, чем их центры тяжести.
Аналитическое уравнение (8.5) показывает, что по мере уменьшения длины стержней R до нуля гравитационный радиус, расстояние между центрами стержней по величине приближается к h. Напротив, увеличение длины стержней с сохранением их массы приводит к увеличению силы их притяжения, поскольку гравитационный радиус уменьшается.
Поведение гравитационных радиусов в общем случае, для произвольных значений длины стержней при неизменной их массе и расстояния между их центрами приведём в виде диаграммы. В блок к аналитическому уравнению (8.5) добавим интегральное уравнение. Используя тот же алгоритм, запишем
После ряда преобразований находим
Теперь для графических построений соединяем оба уравнения в единый блок
В данном варианте мы определяем, чему равно расстояние между точками, гравитационные радиусы, соответствующие разным длинам стержней. Мы просто добавили к расстояниям h индексы: ha - аналитический и hs - интегральный гравитационные радиусы. По смыслу, согласно определению, радиусы для силы между точками FP - это и есть гравитационный радиус. Величины интеграла и логарифмической функции для диаграмм мы уже вычислили для разных значений дискретностей, эквивалентных различным длинам стержней, поэтому теперь мы можем использовать эти же алгоритмы для повторного вычисления. Используя уравнения (8.6), строим графики зависимости гравитационного радиуса от длины стержней
Рис.8.3. Зависимость гравитационного радиуса от длины стержней
На рисунке, на диаграмме вертикальная ось - h, горизонтальная - R. График hg(R) аналитический (синяя штриховая линия) и интегральный (красная линия). Ограниченная точность программного комплекса не позволила построить графики на конечном интервале, вычислить значения уравнений и интеграла вблизи экстремальной точки R = h/2 = 5. Для разрешения этой неопределённости вычислим значения некоторых точек рекурсивно, используя обратный вид уравнения (8.5). Проводим преобразования (8.5) и получаем обратное уравнение в рекурсивном виде
Вид уравнения таков, что непосредственно вычислить радиус R по заданному значению гравитационного радиуса нельзя. Но это уравнение можно решить рекурсией. Задаём некоторое исходное значение радиуса R и находим по этому уравнению его новое значение. Это значение подставляем в это же уравнение и находим следующее значение R.
Вычисляем несколько точек: для недостающего участка графиков рис.8.3 - это точки в трёх нижних строках следующей таблицы. Последнее значение добавлено директивно, поскольку значение конечной точки hg(5) = 0 нам известно.
Две другие точки hg = 8 и hg = 6 мы добавляем для наглядности к полученным выше аналитическим решениям для длин R = h/2, R = h/4 и R → 0. Полученные в таблице решения, значения радиусов R потребовали от 2 до 60 рекурсивных шагов. Наносим эти точки на график рис.8.3. Видим, что аналитические, вычисленные координаты точек совпали с графическими с визуальной, графической точностью. Отметим, что гравитационные центры стержней всегда находятся внутри них. По мере увеличения длины стержней их края сближаются, приближаются к центру системы. Изначально гравитационный центр совпадал с центром стержня, длина которого была близка к нулю. По мере увеличение длины гравитационный центр отстаёт в движении от края стержня, но по мере роста длины стержня, её приближения к половине h, скорость приближения гравитационного центра к центру системы ускоряется, и в последний момент эти центры сливаются. Можно сказать, что гравитационный центр - это условная точка с массой стержня, которая смещается от центра тяжести стержня к центру системы.
Параллельные стержни
Ранее мы рассмотрели задачу о гравитационном притяжении двух параллельных тонких линий. Смысл задачи заключался в сравнении сил притяжения между линиями фиксированной длины на фиксированном удалении друг от друга с силой притяжения между двумя точками той же массы и на том же удалении друг от друга. Теперь мы рассмотрим похожую задачу, но для разного соотношения длины стержней и расстояния между ними. Это соотношение можно менять двояко: либо сближать-удалять два стержня неизменной длины, либо менять длину этих стержней, буквально превращая их в точечные тела. Очевидно, второй вариант математически более удобный.
Рис.8.4. Сравнение силы гравитационного притяжения двух стержней с силой притяжения двух точек такой же массы и на таком же удалении
Массу стержней определим, как и раньше - через их линейную массовую плотность, которая для одинаковых тонких стержней длиной h и массой m определяется уравнением
Понятно, что при изменении длины стержней их линейная плотность также будет меняться. Обозначим стержни символами A и B. Выберем на стержнях элементарные отрезки, "точки" dx и dy. Каждый из этих дифференциальных отрезков, точек имеет массу dm, определяемую произведением длины отрезка на плотность. Соответственно, каждый из этих отрезков одного стержня притягивается к каждому отрезку другого стержня с дифференциальной, малой силой dF
В общем случае эти силы направлены под углом к вертикали. Их проекция на эту вертикаль равна
Здесь переменные x и y - это координаты точек. Полную силу притяжения стержней друг к другу определяем интегрированием
Как и в случае соосных стержней, нас интересует величина силы притяжения в зависимости от длины h стержней. В общем случае величина h может быть любым числом - целым или дробным. Для перехода к целочисленным пределам интегрирования и, соответственно, единичным дифференциалам dx и dy произведём традиционную замену переменных. Диапазоны изменения исходных и переменных для замены соответствуют друг другу.
Связь между новыми и прежними переменными и их дифференциалами имеет вид
Подставляем в уравнение новые переменные и сразу отбрасываем у них штрихи
Выносим дробь h/n за пределы всех скобок и интегралов
Сворачиваем и сокращаем подобные члены
Как обычно мы будем сравнивать силу притяжения стержней с силой притяжения двух точек такой же массы и на таком же удалении друг от друга. Эта сила между точками равна
Вносим изменения в гравитационное уравнение: нижние пределы интеграла меняем на единичные по отмеченным обстоятельствам в VBA-процедуре.
Поскольку мы задались целью сравнения двух сил, обозначим, как и ранее, показатель этого сравнения символом k
Видим, что уравнение содержит безразмерное отношение, которое позволяет упростить запись
Найти табличный интеграл, соответствующий этому уравнению, не удалось. Как мы выявили ранее [3], точность численного решения интеграла зависит от дискретности: чем она выше, тем ближе численное решение к аналитическому решению. Мы сравниваем силу притяжения точек с силой притяжения стержней. Собственно длинам стержней в этом случае мы назначаем длину, кратную их удалённости друг от друга.
Рис.8.5. Сила притяжения двух параллельных стержней
График k(h) на рис.8.5 показывает, что по мере увеличения длины стержней h при неизменном расстоянии d между ними сила их притяжения по сравнению с силой притяжения двух точек с такой же массой стремится по величине к нулю. Уменьшается, стремится к нулю также и абсолютное значение силы притяжения удлиняющихся стержней. Значение асимптоты графика k(0) равно нулю. Напомним, что в рассматриваемых случаях нас интересует не столько величина силы, сколько характер её изменения при изменении длин стержней. При длине стержней, равной нулю, сила их притяжения становится эквивалентна притяжению двух точек той же массы
Это уравнение описывает, в том числе, силу притяжения стержней любой длины, но находящихся на большом удалении друг от друга, то есть, d >> h. Этому условию соответствует график k(d) на рис.8.6, который является коэффициентом, отношением силы притяжения стержней к силе притяжения точек, когда растёт отношение расстояния d между ними к их некоторой фиксированной длине h.
Рис.8.6. Рост относительной силы притяжения двух параллельных стержней при уменьшении их длины
Видим, что этот коэффициент стремится к очевидной единичной асимптоте k(1). Это очевидно, поскольку в пределе такого отношения стержни можно рассматривать как точечные объекты. Это решение показывает, что при расстоянии между стержнями, существенно превышающем их длину, сила притяжения определяется законом Всемирного тяготения Ньютона в его классической формулировке. Добавим, что гравитационные центры стержней во всех рассмотренных случаях параллельных стержней всегда находятся вне области между стержнями.
Таким образом, при любом взаимном расположении тонких стержней сила их гравитационного притяжения описывается классическим законом Всемирного тяготения только в случае, когда их длина много меньше расстояния между ними.
9. График параметрической интегральной функции
Как мы выше определили, постоянный интеграл - это число. Именно число, а конкретно в рассмотренных в данной работе случаях - коэффициент относительной силы, который был главной целью вычисления основной массы интегралов. Достоверность, корректность этого числа, коэффициента мы определяли по тенденции интеграла и предельному значению дискретности, её минимальному значению, при котором интеграл сходился. Но постоянный параметрический интеграл - это уже массив чисел, которые можно представить как в виде таблицы, так и в виде графика. Параметры интегралов могут входить как непосредственно в его подынтегральную функцию, так и в пределы интегрирования. Самым простым параметрическим интегралом является, очевидно, однократный интеграл.
В таком виде мы рассмотрели функцию с интегралом косинуса в примере на рис.1.1 и в примере по ссылке рис.1.4. Хотя однократный интеграл даёт площадь под интегрируемой функцией, параметрические пределы интегрирования уже можно рассматривать как отдельные точки такого вычисленного интеграла. Более интересным и сложным является пример кратных интегралов. В случае их параметрического представления, особенно, в случае параметрических пределов интегрирования, эти интегралы дают нам объект в многомерном пространстве соответствующей размерности. Явную геометрическую функцию, функцию в виде графика даёт нам интеграл однократный параметрический. Производная от этой функции является подынтегральной функцией. То есть, такой интеграл - это не одно число, а двухмерная таблица или плоский график. Двойной интеграл - это уже трёхмерная таблица или некая поверхность. Рассмотрим ещё более сложный вариант - четырёхкратный интеграл (5.4), сформированный в задаче о притяжении двух зеркал.
Полное параметрическое представление его пределов интегрирования, очевидно, даст на функцию четырёх переменных F(u, v, y, x), графическое изображение которой, видимо, представляет собой объёмное массивное тело переменной плотности. Действительно, параметрический одинарный интеграл графически даёт нам линию. Двойной - поверхность. Тройной - объёмное тело. Следует предположить, что четверной интеграл придаёт каждой из точек объёмного тела некую массу. Впрочем, существует мнение, что такой интеграл даёт объёмный четырёхмерный объект типа гиперкуба (тессеракта).
Установим в интеграле (9.1) параметрические пределы интегрирования. Для этого рассмотрим на верхнем зеркале небольшой прямоугольный участок в его центре, а на нижнем - линию небольшой ширины вдоль оси x. Относительная ортогональная сила притяжения этого участка к линии, коэффициент k можно описать параметрическими пределами интегрирования
Координаты, границы этих прямоугольных фигур удобно задать через соответствующие значения дискретности каждой из переменных.
В этом наборе параметров действительно переменной является только параметр ny. Изменение этого параметра означает "перемещение" полосы на нижнем зеркале. Сам интеграл (9.1), следовательно, даёт нам относительные ортогональные силы притяжения участка на верхнем зеркале к разным полоскам на нижнем. Напомним, что относительная сила притяжения - это коэффициент k, показывающий, насколько, во сколько раз сила притяжения зеркал или его участков отличается от силы притяжения двух точек той же массы и на том же удалении друг от друга. Ортогональность означает составляющую силы, перпендикулярную к плоскости зеркал. Без нарушения общности установим Δu = Δv = Δy = Δ. Тогда можно сказать, что параметры u0u1 - это ширина площадки на верхнем зеркале, v0v1 - его длина. Параметры y0y1 - это ширина полоски на нижнем зеркале, а параметр ny - её координата. Соответственно, параметр x задаёт полную длину полоски от 1 до n. Для формирования компьютерной процедуры заменим символ Δ печатным символом m. Устанавливаем константы n = 1'000; Δ = m = 5. С этими обозначениями интеграл в общем виде можно записать следующим образом
Подставляем значения условных переменных, пределов в интеграл, который окончательно приобретает вид
Используя это уравнение, строим график kuv(ny).
Рис.9.1. Сила притяжения площадки на верхнем зеркале к полоске на нижнем в зависимости от её положения
Строго говоря, картину можно назвать очевидной. Наибольшая сила притяжения площадки к полоске возникает в том месте, в котором эта площадка на верхнем зеркале находится точно над полоской нижнего зеркала. Притяжение площади к полоскам, находящимся на краях нижнего зеркала, минимально.
Выводы
Численное интегрирование позволяет произвести вычисление определённых интегралов, не являющихся табличными, либо требующих сложных преобразований для сведения к табличным. Если интеграл, подынтегральная функция является параметрической, то есть, функцией от некоторой переменной, то численное интегрирование позволяет построить интегральный график как функцию этой переменной.
Численное интегрирование применимо к интегралам любой кратности: двойным, тройным и выше. Однако при этом время вычисления интеграла может быть большим. В некоторых случаях время можно уменьшить интервальным интегрированием, то есть, вычислением интеграла на части области определения аргументов и последующим пропорциональным, кратным расширением решения на всю область определения.
Точность вычисления определяется дискретностью переменных подынтегральной функции. Чем выше дискретность, тем точнее результат. Существует некоторая предельная дискретность, превышение которой не приводит к существенному повышению точности, то есть, является границей сходимости интеграла.
Гравитационное притяжение объёмных тел в общем случае не описывается точно законом всемерного тяготения Ньютона в классической формулировке. Расчёты с использованием численного интегрирования показали, что точка притягивается не к центру объёмного тела, а к точке, находящейся чуть дальше его центра. Чем дальше точка от такого тела, тем ближе к его центру эта эквивалентная точка притяжения, которую, в отличие от центра тяжести или геометрического центра, можно назвать гравитационным центром. Положение гравитационного центра зависит от положения притягиваемого тела и для разных объектов находится в разных местах.
Исключением являются сферические тела - однородные шары, центром притяжения, гравитационным центром которых является их центр тяжести, геометрический центр. Притяжение шаров друг к другу или точки к шару законом Всемирного тяготения Ньютона в классической форме описывается точно.
Сила гравитационного притяжения двух плоских тел меньше силы притяжения двух точек, имеющих ту же массу и находящихся на таком же удалении друг от друга. Такое же соотношение наблюдается и для массивных линий, отрезков.
ПРИЛОЖЕНИЯ
П1. Excel - процедура
Приведён сокращённый текст компьютерной процедуры, использованной при численном интегрировании уравнения (5.10). Это основной блок процедуры без вспомогательных процедур и деклараций.
Sub numbinteg()
'Базовые параметры интегралов
k_masht = Sheets("track").Range("mashtab")
positie = Sheets("track").Range("position")
'Обнуляем для установки первого значения времени шага
t_step3 = 0
b_step = 0
period_info = Sheets("track").Range("period_info")
'Установка размеров и положения диаграммы
diagr = 1115
d_Top = 0
d_Width = 540
d_Left = diagr - d_Width
d_Height = 340
'Имя диаграммы внести вручную
With Worksheets("track").Shapes("Диаграмма 11")
.Left = d_Left
.Top = d_Top
.Width = d_Width
.Height = d_Height
End With
'Установим заголовки таблицы
Sheets("track").Range("tendenc") = "Тенднц."
Sheets("track").Range("integral") = "Интеграл"
'Эти три ячейки должны быть в этом порядке
Sheets("track").Range("time_proshlo") = "Прошло"
Sheets("track").Range("time_end") = "Завершение"
Sheets("track").Range("step_time") = "Шаг"
Sheets("track").Range("predel") = "n%"
'Заполнение колонки selekt только при sel_start = 0
sel_start = Sheets("track").Range("selekt")
'Дискретности начинаем считывать с ячейки sel = 1
'Процедура ТОЛЬКО для первичной установки параметров
'При значении sel_start = 0 процедура numbinteg не запускается
If sel_start = 0 Then
sel_start = 1
Sheets("track").Range("selekt") = sel_start
'Обнуляем для установки первого значения времени шага
t_step3 = 0
b_step = 0
Sheets("track").Range("mods").Offset(2) = 0
b_step_y = 0
y_mirror = 0
Sheets("track").Range("mods") = ""
'Первичное значение шага в таблице Скачок
step_imp = 1
Sheets("track").Range("mods").Offset(3) = 1
'Заполняем колонку selekt значениями v_start
v_table = 2
For i = 1 To 14
v_table = v_table * 2
Sheets("track").Range("selekt").Offset(i) = v_table
Next
'Последнее значение для аварийного выхода - пустая ячейка
Sheets("track").Range("selekt").Offset(i) = ""
'Выходим, поскольку блок только для формирования списка
Exit Sub
End If
'Заголовки таблицы Скачка
Sheets("track").Range("impact") = "Скачок"
Sheets("track").Range("impact").Offset(0, 1) = " Тенденция_" & positie & "%"
Tk_index = Sheets("track").Range("selekt").Offset(sel_start)
Sheets("track").Range("impact").Offset(0, 2) = " Интеграл_" & Tk_index
'ГЛАВНЫЙ ЦИКЛ перебора ДИСКРЕТНОСТЕЙ
'Значение sel равно номеру ячейки с дискретностью, начиная с 1
'Первый запуск происходит со значением sel = 1
For sel = sel_start To 100
'Считываем текущее значение v_end = 100n / positie
v_end = Sheets("track").Range("selekt").Offset(sel)
'Если считано Empty, конец таблицы, выходим из цикла
'Это точка реального завершения процедуры
If v_end = "" Then
Call Beep_melody
'Напоминание о необходимости сохранения файла
MsgBox ("Вычисление интеграла закончено." & Chr(13) & _
"Сохранить файл - обязательно!")
'Вычисление завершено, удаляем окраску ячейки
Sheets("track").Range("selekt").Offset(sel, 1). _
Interior.Pattern = xlNone
'Это основная точка выхода, завершения процедуры
Exit Sub
End If
'Устанавливаем окраску текущей ячейки
Sheets("track").Range("selekt").Offset(sel, 1). _
Interior.color = &HFFFF&
If sel >= 2 Then
'Снимаем окраску c предыдущей строки. Заодно и со следующей
Sheets("track").Range("selekt").Offset(sel - 1, 1). _
Interior.Pattern = xlNone
Sheets("track").Range("selekt").Offset(sel + 1, 1). _
Interior.Pattern = xlNone
End If
'Сохраняем значение n
n = v_end * 100
'При установленном в задаче интервале = 1 доля от n
v_end = n * positie / 100
v_start = v_end
'Проверяем значение главного счётчика
'Перед пуском он обязательно должен быть задан
mods_start = Sheets("track").Range("mods")
'Если он пустой, Empty, или = 0, то это первый старт
If mods_start = "" Or mods_start = 0 Then
'Задаём исходные, стартовые значения, параметров
'Очистка области диаграммы Скачка
ActiveWorkbook.Names.Add Name:="impact_start", RefersTo:=Range("impact").Offset(1)
ActiveWorkbook.Names.Add Name:="impact_end", RefersTo:=Range("impact").Offset(8000, 3)
Range("impact_start:impact_end").Select
Selection.ClearContents
ku_mirror = 0
Sheets("track").Range("mods").Offset(2) = 0
Sheets("track").Range("integral").Offset(sel) = 1
Sheets("track").Range("tendenc").Offset(sel) = "расчет..."
Sheets("track").Range("time_proshlo").Offset(sel) = ""
Sheets("track").Range("time_end").Offset(sel) = ""
Sheets("track").Range("step_time").Offset(sel) = ""
step_imp = 1
Sheets("track").Range("mods").Offset(3) = step_imp
Sheets("track").Range("mods") = 0
End If
'Переводим счётчик шагов, значение mods на следующий шаг
mods_start = Sheets("track").Range("mods") + 1
Sheets("track").Range("mods") = mods_start
ku_mirror = Sheets("track").Range("mods").Offset(2)
'Считываем время старта.
'Предварительно нужна запись чего угодно куда угодно,
'иначе из (0, -1) считывалось прошлое значение
time_start = Sheets("track").Range("mods").Offset(0, -1)
Sheets("track").Range("mods").Offset(4) = time_start
step_proshlo = 0
step_time = time_start
'Смещаем экран в начальное положение
Sheets("track").Range("mods").Offset(0, -1).Activate
Sheets("track").Range("integral").Offset(sel).Activate
Columns("a:w").EntireColumn.AutoFit
'ШШШШШШ Н А Ч А Л О процедуры ШШШШШШ
step_imp = Sheets("track").Range("mods").Offset(3)
p = d_mirror / h_mirror
'Групповой параметр: добавлен множитель n
p3n = p ^ 3 * n
'Параметры для вычисления среднего времени t_process
n_proc_start = b_step
t_proc_start = 0
'ЦИКЛ С_Т_А_Р_Т С_Т_А_Р_Т С_Т_А_Р_Т С_Т_А_Р_Т
'Интегрирование от v_start до v_end = v_start
'================== v_start
For v_mirror = v_start To v_end
'Начальное значение счетчика внутреннего цикла по y
b_step_y = 0
'Продолжение вычислений, если был останов
y_mirror_step = Sheets("track").Range("mods")
'================== y_start
For y_mirror = y_mirror_step To n
'================== x_start
vyp = (v_mirror - y_mirror) ^ 2 + n ^ 2 * p ^ 2
For x_mirror = 1 To n
k2 = 1 / ((1 - x_mirror) ^ 2 + vyp) ^ (3 / 2)
ku_mirror = ku_mirror + k2
Next 'x
'================== x_end
'Счетчик по y
b_step_y = b_step_y + 1
'Счётчик шагов для измерения времени шага
step_proshlo = step_proshlo + 1
If b_step_y >= period_info Then
b_step_y = 0
Sheets("track").Range("mods") = y_mirror
'Время выводим с каждым выводом y_mirror
time_stop = Sheets("track").Range("mods").Offset(0, -1)
Application.StatusBar = "Шаг = " & y_mirror
'Активируем ячейку Интеграл
Sheets("track").Range("integral").Offset(sel).Activate
'Время выводим с каждым выводом y_mirror
time_proshlo = time_stop - time_start
Sheets("track").Range("time_proshlo").Offset(sel) = time_proshlo
'Определяем среднее время шага за период period_info
'После последнего измерения прошло время
step_time = time_stop - step_time
'Вносим в таблицу найденное время шага
Sheets("track").Range("step_time").Offset(sel) = step_time
'Осталось шагов до окончания после перезапуска
step_to_end = n - y_mirror
'Время до завершения текущее, остаток
time_ostatok = step_time * step_to_end / period_info
Sheets("track").Range("time_end").Offset(sel) _
= time_stop + time_ostatok
Sheets("track").Range("mods").Offset(5) = time_ostatok
'Запоминаем новое время начала отсчёта шагов
step_time = time_stop
'Заносим в ячейки Интеграл все промежуточные значения
kx_end_2 = ku_mirror * p3n
ku_mirror_2 = kx_end_2 * k_masht
Sheets("track").Range("integral").Offset(sel) = ku_mirror_2
'Сохраняем текущие значения для продолжения после останова
Sheets("track").Range("mods").Offset(2) = ku_mirror
Sheets("track").Range("mods") = y_mirror
'Устанавливаем значение Ts - тенденции
integ = Sheets("track").Range("integral").Offset(sel)
integ2 = Sheets("track").Range("integral").Offset(sel - 1)
'При отсутствии предыдущего, самого верхнего значения
If sel = 1 Then integ2 = integ * 4
tend = integ2 / integ
Sheets("track").Range("tendenc").Offset(sel) = tend
'Заполняем таблицу Скачка. Шаг таблицы в %
imp_percent = 100 * y_mirror / n
Sheets("track").Range("impact").Offset(step_imp) = imp_percent
Sheets("track").Range("impact").Offset(step_imp, 1) = tend
Sheets("track").Range("impact").Offset(step_imp, 2) = ku_mirror_2
step_imp = step_imp + 1
Sheets("track").Range("mods").Offset(3) = step_imp
Columns("a:w").EntireColumn.AutoFit
End If
Next 'y
'================== y_end
b_step = b_step + 1
kx_end_2 = ku_mirror * p3n
ku_mirror_2 = kx_end_2 * k_masht
Sheets("track").Range("integral").Offset(sel) = ku_mirror_2
'При отладке Activate могла сместиться. Восстанавливаем
Application.StatusBar = "Шаг = " & y_mirror
Sheets("track").Range("integral").Offset(sel).Activate
Columns("a:w").EntireColumn.AutoFit
Next 'v
'================== v_end
'ЦИКЛ К_О_Н_Е_Ц ЦИКЛ. ЦИКЛ. ЦИКЛ. ЦИКЛ.
Call Beep_melody3
'Выводим окончательное значение интеграла
'Цикл может закончиться без захода в эту процедуру
kx_end_2 = ku_mirror * p3n
ku_mirror_2 = kx_end_2 * k_masht
Sheets("track").Range("integral").Offset(sel) = ku_mirror_2
'ШШШ К О Н Е Ц П Р О Ц Е Д У Р Ы numbinteg ШШ
'Переносим в таблицу дату и время завершения строки
date_end = Sheets("track").Range("mods").Offset(0, -1)
Sheets("track").Range("time_end").Offset(sel) = date_end
'Вычисляем время, затраченное на расчёт строки
time_proshlo = date_end - time_start
'Записываем затраченное время в расчётную ячейку
Sheets("track").Range("time_proshlo").Offset(sel) = time_proshlo
'Устанавливаем окончательно значение Ts - тенденции
integ = Sheets("track").Range("integral").Offset(sel)
integ2 = Sheets("track").Range("integral").Offset(sel - 1)
'При отсутствии предыдущего, самого верхнего значения
If sel = 1 Then integ2 = integ * 4
tend = integ2 / integ
Sheets("track").Range("tendenc").Offset(sel) = tend
'Увеличиваем на 1 номер ячейки для расчёта - следующей ячейки
sel_start = sel_start + 1
'Сохраняем это значение в ячейке selekt
Sheets("track").Range("selekt") = sel_start
Columns("a:w").EntireColumn.AutoFit
Next 'sel
'До этого места процедура не дойдёт никогда.
'Приведена для законченности
'Все дискретности из таблицы РАССЧИТАНЫ
Call Beep_melody
'Напоминание о необходимости сохранения файла
MsgBox ("Вычисление интеграла закончено." & Chr(13) & "Сохранить файл - обязательно!")
End Sub
Хотя процедура и выглядит относительно объёмной, на самом деле её ядром, главным "вычислительным центром" являются всего лишь следующие несколько строк:
For v_mirror = v_start To v_end
For y_mirror = y_mirror_step To n
vyp = (v_mirror - y_mirror) ^ 2 + n ^ 2 * p ^ 2
For x_mirror = 1 To n
k2 = 1 / ((1 - x_mirror) ^ 2 + vyp) ^ (3 / 2)
ku_mirror = ku_mirror + k2
Next 'x
Next 'y
Next 'v
Остальные строки - вспомогательные, информационные, комментарии. На следующем рисунке представлен скриншот рабочего стола при вычислении интеграла (5.10)
Рис.П1. Скриншот рабочего стола
Приведённый текст процедуры и иллюстрация предназначены показать, что рассмотренная задача является рядовой компьютерной процедурой и не содержит ничего особого сложного.
П2. Загадочное уравнение
В процессе вычисления уравнений в задаче о притяжении точки к диску обнаружено странное уравнение. Решение его логически известно: коэффициент, интеграл аналитический и числовой при росте параметра от нуля до бесконечности должен стремиться к единице. Аналитическое решение имеет компактный вид, выглядит как обычное уравнение
Кажется, что решение его не представляет особого труда. Действительно, проделаем несколько тривиальных преобразований. Поскольку они достаточно краткие, приведём их полностью, чтобы возможные некорректные преобразования были предельно открыты
И тут возникает странная ситуация. Проверка этого решения, полученного вполне корректными преобразованиями, не подтверждает корректность этого решения. Действительно, подставляем его в уравнение, из которого оно получено, и делаем все необходимые преобразования
Видим, что найденное решение уравнения его решением не является. Вместе с тем, точное решение уравнения нам известно: это величина параметра p, стремящаяся к бесконечности. Как вариант, величина p равна очень большому числу, но никак не полученному в процессе рассмотренного решения.
Ссылки
1. Путенихин П.В. Логические основания многомерных пространств. -- Саратов: "АМИРИТ", 2018. - 396 с., цв. илл., URL:
https://www.elibrary.ru/item.asp?id=42690781
2. Путенихин П.В. Численное интегрирование, 2022, URL:
http://samlib.ru/p/putenihin_p_w/numinteg.shtml
http://samlib.ru/img/p/putenihin_p_w/numinteg/numinteg455.pdf
3. Дифференциал площади круга dS, 2021, URL
http://samlib.ru/editors/p/putenihin_p_w/ds_199.shtml
04.28 - 14.08.2022