CFA LogoCFA Logo Computer
Новости Статьи Магазин Прайс-лист Драйвера Контакты
Новости
RSS канал новостей
В конце марта компания ASRock анонсировала фирменную линейку графических ускорителей Phantom Gaming. ...
Компания Huawei продолжает заниматься расширением фирменной линейки смартфонов Y Series. Очередное ...
Компания Antec в своем очередном пресс-релизе анонсировала поставки фирменной серии блоков питания ...
Компания Thermalright отчиталась о готовности нового высокопроизводительного процессорного кулера ...
Компания Biostar сообщает в официальном пресс-релизе о готовности флагманской материнской платы ...
Самое интересное
Программаторы 25 SPI FLASH Адаптеры Optibay HDD Caddy Драйвера nVidia GeForce Драйвера AMD Radeon HD Игры на DVD Сравнение видеокарт Сравнение процессоров

АРХИВ СТАТЕЙ ЖУРНАЛА «МОЙ КОМПЬЮТЕР» ЗА 2003 ГОД

Мысли о Паскале

Владислав ДЕМЬЯНИШИН nitromanit@mail.ru

Продолжение, начало см. в МК №46, 51-52, 4, 6-7, 10, 12-13, 16-18, 22, 24, 29, 34, 41, 46, 4, 6, 17, 21, 23, 28, 30, 32, 39 (165, 170-171, 175, 177-178, 181, 183-184, 187-189, 193, 195, 200, 205, 212, 217, 227, 229, 240, 244, 246, 251, 253, 255, 262).

Спрашивали? Отвечаю…

Расширение математических возможностей Паскаля

Однажды, когда мне понадобилось реализовать сложный алгоритм с использованием логарифмов, возведений в произвольную степень и вычислением определенных интегралов, я был неприятно удивлен, что среди встроенных алгебраических функций среды Turbo Pascal подобные возможности отсутствуют. Тогда, поворошив в памяти остатки еще не выветрившихся знаний из школьного и вузовского курсов алгебры и высшей математики, я набросал исходный код нескольких необходимых функций. С подбором эффективного алгоритма вычисления определенного интеграла мне помогла книга, указанная в списке литературы в конце этой статьи. Итак, сегодня предлагаю вам составить модуль math.pas, который мы наделим солидными математическими возможностями.

Как обычно, начнем с заголовка модуля и опишем его традиционно, описав в интерфейсной части тип, который нам понадобится для реализации обратного вызова из функции вычисления определенного интеграла TIntegralFunc. Надо сказать, что у меня в модуле везде использовался тип Real, но когда я сел за написание данной статьи, мне захотелось сделать модуль универсальным применительно к различным вещественным типам и заодно избавить вас от тупой работы — сами понимаете, насколько утомительно рыскать по коду, повсюду меняя базовый тип модуля на Single или Double, по мере необходимости. Поэтому давайте опишем тип Float, который, если что, легко можно будет образовать от любого другого вещественного типа, исправив лишь одну строчку в описании type. Ну и какая же интерфейсная часть модуля обойдется без перечня заголовков публикуемых подпрограмм?

А теперь настал черед раздела implementation, в котором сформируем исходный код для каждой подпрограммы. Начнем с простых, но в то же время наиболее востребованных.

Логарифмы

Так как разработчики Turbo Pascal дали нам очень скромный джентльменский набор функций Ln и Exp для вычисления натурального логарифма и экспоненты, то для нахождения логарифма числа X по основанию a (то бишь Log a X) воспользуемся одним из свойств этого логарифма — последний будет равен отношению логарифма числа X по некоторому основанию b к логарифму числа a по некоторому основанию b; главное, чтобы в обоих логарифмах дроби основание было одинаковое:

Таким образом, в качестве основания логарифмов дроби можно взять экспоненту:

и тогда (по свойству логарифма Log e X = Ln X) составим иную формулу вычисления логарифма числа X по основанию a:

В итоге получаем функцию:

Это же свойство сгодится и для нахождения десятичного логарифма числа X (то бишь Lg X по основанию 10):

Чтобы уменьшить количество вычислений в данной формуле, предварительно найдем значение константного выражения 1/Ln(10) и занесем его в переменную в основном блоке Begin..End модуля:

теперь формула будет как минимум на одно действие проще:

В итоге получим функцию

которая позволит вычислять десятичный логарифм от X.

Степени

Теперь поставим перед собой задачу возвести число X в степень N (для функции SqrN), для чего воспользуемся следующим свойством логарифма (пускай знак ^ указывает на возведение в степень):

Ничто не мешает нам в качестве основания логарифма a выбрать экспоненту (e)

тогда по свойству логарифма Log e X = Ln X данное равенство можно представить так:

Теперь необходимо вспомнить, как звучит определение логарифма. Логарифмом числа X по основанию a называется такой показатель степени b, в который надо возвести основание a, чтобы получить число X. То есть, имеем свойства логарифма Log a X = b и a^b = X, из которых следует формула:

Тогда функция SqrN должна выглядеть так:

Последняя задача, которая может касаться логарифмов — это извлечение корня N-степени из числа X (для функции SqrtN). Для решения этой задачи следует вспомнить свойства степеней, из которых следует, что корень N-степени из числа X равен возведению числа X в степень 1/N. Тогда, используя формулу предыдущей функции, получим:

Здесь мне не нравится деление 1/N — его можно сократить, пользуясь свойством дробей, и тогда окончательный вид составляемой функции SqrtN будет таков:

На всякий случай определение знака функцией Sign, которая в случае положительного аргумента возвращает 1, а в случае отрицательного —–1.

Тригонометрия

Ну и немножко тригонометрии, функции которой, думаю, не нуждаются в комментариях. Что? Нуждаются? Ладно, мне это лишь доставит удовольствие.

Я не сделаю революционного открытия Америки через форточку, если скажу, что в тригонометрии арксинус вычисляется по формуле:

arcsin = arctg (X/SQRT(1–X*X))

Формула довольно проста, но при ее использовании можно наткнуться на «подводные камни» — если в качестве аргумента X будет задана величина 1, то выражение SQRT(1–X*X) даст 0, а так как на ноль делить нельзя, то справедливо получим замечание в виде соответствующей ошибки исполнения. Чтобы избежать такого казуса, я решил прибегнуть к системе упреждения ошибки, предварительно вычислив выражение 1–X*X — если оно дает нулевой результат, то функция возвращает результат PI/2 со знаком аргумента.

Для вычисления арккосинуса можно применить формулу:

arccos = arctg (SQRT(1–X*X)/X)

Но и она не лишена коварных «коралловых рифов», так как при нулевом значении аргумента моментально приведет к ошибке деления на ноль. Поэтому приходится предварительно проверять значение аргумента, и если оно равно нулю, то функция возвращает результат PI/2.

Итак, жаждете получить арксинус? Их есть у меня :-)!

Вызывайте функцию ArcSin:

и ее сестричка, функция ArcCos:

И, конечно, их внучатая племянница, функция ArcCtg:

Вычисление определенного интеграла

Вот теперь рассмотрим алгоритм вычисления определенного интеграла. Основная задача численного интегрирования сводится к нахождению значения собственного определенного интеграла, подынтегральная функция которого на отрезке [a, b] не имеет особенностей.

В общем случае интервал интегрирования [a, b] разбивается на M частей. В свою очередь каждая из них делится на N частей, в пределах каждой из которой y=f(x) аппроксимируется полиномом, интегрирование которого возможно по достаточно простым формулам.

Я выбрал алгоритм численного интегрирования методом парабол (Симпсона), который относится к числу простых методов интегрирования, но дает наиболее высокую точность, потому и чаще всего применяется. Данный метод позволяет употребить переменный шаг, выбираемый автоматически из условия получения заданной точности E результата. Для этого величине P2 придаются значения 2, 4, 8, 16 и т.д. При каждом удвоении P2 точность улучшается приблизительно в 15 раз. Процесс интегрирования следует прекратить при выполнении условия (Fi(X)–Fi+1(X))/15<E.

При подготовке к интегрированию вычисляется коэффициент точности T=Sqr(15*E), где роль коэффициента E выполняет параметр Accuracy. Следует учесть, что максимальная точность ограничена размером мантиссы вещественного типа, от которого образован тип Float, то есть при коэффициенте точности, например, равном 0.000000000000000000001 (т.е. 1.0E–21 или, проще говоря, до 21 знака после запятой) и типе Real выражение T:=Sqr(15*Accuracy) даст нулевой результат, что будет значить переполнение данного типа, и это приведет к зацикливанию алгоритма интегрирования. Поэтому следует либо определять точность до конкретного разряда после запятой, либо указать 0 в качестве параметра Accuracy, и тогда функция вычисления интеграла сама установит максимально допустимую точность для типа, от которого образован тип Float. В коде функции также выполняется защита от отрицательной величины Accuracy. Параметры A и B задают интервал интегрирования подынтегральной функции Fx, которой может быть функция с дальним вызовом Far, имеющая заголовок function( X : float ) : float:

Рассмотрим применение данной функции на примере:

Легко убедиться, что теперь вычисление определенного интеграла превратится в простую и наглядную операцию.

Чем выше заданная точность (т.е. меньше значение параметра Accuracy), тем больше время интегрирования. В большинстве случаев задаваемую точность можно свести к минимуму, увеличив значение параметра Accuracy, что позволит получить весьма ощутимый прирост в производительности при вычислении интеграла.

Чтобы сократить время вычисления при достаточно высокой точности, следует использовать более сложные методы численного интегрирования, например методы Бодэ, Ньютона-Котеса, Уэддля, Чебышева, Гаусса и др.

Пришел, увидел, а побеждать уже нечего. :-).

(Продолжение следует)

Литература

1. Дьяконов В.П. Справочник по расчетам на микрокалькуляторах – М.:Наука, 1989. – 462 с.

Рекомендуем ещё прочитать:






Данную страницу никто не комментировал. Вы можете стать первым.

Ваше имя:
Ваша почта:

RSS
Комментарий:
Введите символы или вычислите пример: *
captcha
Обновить





Хостинг на серверах в Украине, США и Германии. © www.sector.biz.ua 2006-2015 design by Vadim Popov