Продолжение, начало см. в МК №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. Ну и какая же интерфейсная часть модуля обойдется без перечня заголовков публикуемых подпрограмм?
Unit Math; interface type Float = Real; TIntegralFunc = function( X : float ) : float; var Ln10 : float; function Integral( A, B, Accuracy : float; Fx : TIntegralFunc ) : float; function Lg( X : float ) : float; function Log( A, N : float) : float; function SqrN( X, N : float ) : float; function SqrtN( X, N : float ) : float; function ArcSin( X : float ) : float; function ArcCos( X : float ) : float; function ArcCtg( X : float ) : float;
А теперь настал черед раздела implementation, в котором сформируем исходный код для каждой подпрограммы. Начнем с простых, но в то же время наиболее востребованных.
Логарифмы
Так как разработчики Turbo Pascal дали нам очень скромный джентльменский набор функций Ln и Exp для вычисления натурального логарифма и экспоненты, то для нахождения логарифма числа X по основанию a (то бишь Log a X) воспользуемся одним из свойств этого логарифма — последний будет равен отношению логарифма числа X по некоторому основанию b к логарифму числа a по некоторому основанию b; главное, чтобы в обоих логарифмах дроби основание было одинаковое:
Log a X = Log b X / Log b a;
Таким образом, в качестве основания логарифмов дроби можно взять экспоненту:
Log a X = Log e X / Log e a;
и тогда (по свойству логарифма Log e X = Ln X) составим иную формулу вычисления логарифма числа X по основанию a:
Log a X = Ln X / Ln a;
В итоге получаем функцию:
function Log( a, X : float) : float; begin Log := Ln(X) / Ln(a); end;
Это же свойство сгодится и для нахождения десятичного логарифма числа X (то бишь Lg X по основанию 10):
Lg X = Log 10 X = Ln X / Ln 10;
Чтобы уменьшить количество вычислений в данной формуле, предварительно найдем значение константного выражения 1/Ln(10) и занесем его в переменную в основном блоке Begin..End модуля:
begin Ln10 := 1/Ln(10); end.
теперь формула будет как минимум на одно действие проще:
Lg X = Ln10 * Ln X;
В итоге получим функцию
function Lg( X : float ) : float; begin Lg := Ln10*Ln(X); end;
которая позволит вычислять десятичный логарифм от X.
Степени
Теперь поставим перед собой задачу возвести число X в степень N (для функции SqrN), для чего воспользуемся следующим свойством логарифма (пускай знак ^ указывает на возведение в степень):
Log a X^N = N * Log a X;
Ничто не мешает нам в качестве основания логарифма a выбрать экспоненту (e)
Log e X^N = N * Log e X;
тогда по свойству логарифма Log e X = Ln X данное равенство можно представить так:
Ln X^N = N * Ln X;
Теперь необходимо вспомнить, как звучит определение логарифма. Логарифмом числа X по основанию a называется такой показатель степени b, в который надо возвести основание a, чтобы получить число X. То есть, имеем свойства логарифма Log a X = b и a^b = X, из которых следует формула:
X^N = e ^ ( N * Ln X );
Тогда функция SqrN должна выглядеть так:
function SqrN( X, n : float ) : float; begin SqrN := Exp( n*Ln(X) ); end;
Последняя задача, которая может касаться логарифмов — это извлечение корня N-степени из числа X (для функции SqrtN). Для решения этой задачи следует вспомнить свойства степеней, из которых следует, что корень N-степени из числа X равен возведению числа X в степень 1/N. Тогда, используя формулу предыдущей функции, получим:
SqrtN = Exp( 1 / N * Ln(X) );
Здесь мне не нравится деление 1/N — его можно сократить, пользуясь свойством дробей, и тогда окончательный вид составляемой функции SqrtN будет таков:
function SqrtN( X, n : float ) : float; begin SqrtN := Exp( Ln(X) / n ); end;
На всякий случай определение знака функцией Sign, которая в случае положительного аргумента возвращает 1, а в случае отрицательного —–1.
function Sign( X : float ): float; begin if X = 0 then Sign := 1 else Sign := X / abs(X); end;
Тригонометрия
Ну и немножко тригонометрии, функции которой, думаю, не нуждаются в комментариях. Что? Нуждаются? Ладно, мне это лишь доставит удовольствие.
Я не сделаю революционного открытия Америки через форточку, если скажу, что в тригонометрии арксинус вычисляется по формуле:
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:
function ArcSin( X : float ) : float; var y, s : float; begin y := 1 – X * X; y := abs( y ); s := sign( X ); if y = 0 then ArcSin := s * PI / 2 else ArcSin := ArcTan( X / Sqrt( y ) ); end;
и ее сестричка, функция ArcCos:
function ArcCos( X : float ) : float; var y, s : float; begin y := 1 – X * X; y := abs( y ); s := sign( X ); if X = 0 then ArcCos := PI / 2 else ArcCos := ArcTan ( Sqrt( y ) / X ); end;
И, конечно, их внучатая племянница, функция ArcCtg:
function ArcCtg( X : float ) : float; begin ArcCtg := PI / 2 – ArcTan( X ); end;
Вычисление определенного интеграла
Вот теперь рассмотрим алгоритм вычисления определенного интеграла. Основная задача численного интегрирования сводится к нахождению значения собственного определенного интеграла, подынтегральная функция которого на отрезке [a, b] не имеет особенностей.
В общем случае интервал интегрирования [a, b] разбивается на M частей. В свою очередь каждая из них делится на N частей, в пределах каждой из которой y=f(x) аппроксимируется полиномом, интегрирование которого возможно по достаточно простым формулам.
Я выбрал алгоритм численного интегрирования методом парабол (Симпсона), который относится к числу простых методов интегрирования, но дает наиболее высокую точность, потому и чаще всего применяется. Данный метод позволяет употребить переменный шаг, выбираемый автоматически из условия получения заданной точности E результата. Для этого величине P2 придаются значения 2, 4, 8, 16 и т.д. При каждом удвоении P2 точность улучшается приблизительно в 15 раз. Процесс интегрирования следует прекратить при выполнении условия (Fi(X)–Fi+1(X))/15= 0; p1 := p1 + C; Fold := F; F := ( (p1 + C) / 3 )*BA; until Sqr( Fold – F ) –T < 0; Integral := F; end;`
Рассмотрим применение данной функции на примере:
Uses Math; function ifunc( x : float ) : float; far; begin { где f(x[0..1]) = 1.3987174742 } ifunc := Sqrt( 2*x + 1 ); end; function ifunc2( x : float ) : float; far; begin { где f(x[1..5]) = 0.9074539 } ifunc2 := (x*x*x) / (x*x*x*x + 16); end; begin writeln('Calc ',Integral(0,1,0,ifunc):1:10); writeln('Calc2 ',Integral(1,5,0.01,ifunc2):1:10); end.
Легко убедиться, что теперь вычисление определенного интеграла превратится в простую и наглядную операцию.
Чем выше заданная точность (т.е. меньше значение параметра Accuracy), тем больше время интегрирования. В большинстве случаев задаваемую точность можно свести к минимуму, увеличив значение параметра Accuracy, что позволит получить весьма ощутимый прирост в производительности при вычислении интеграла.
Чтобы сократить время вычисления при достаточно высокой точности, следует использовать более сложные методы численного интегрирования, например методы Бодэ, Ньютона-Котеса, Уэддля, Чебышева, Гаусса и др.
Пришел, увидел, а побеждать уже нечего. :-).
(Продолжение следует)
Литература
1. Дьяконов В.П. Справочник по расчетам на микрокалькуляторах – М.:Наука, 1989. – 462 с.