Автор: Владислав Демьянишин

 

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

 

Так что сегодня предлагаю вам составить модуль MATH.PAS, который мы наделим, прямо скажем, нехилыми вычислительными возможностями в области математики.

 

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

  

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

 

Unit Math;

 

 

       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 по основанию a, который равен отношению логарифма числа 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;

Log := Ln(X) / Ln(a);

 

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

 

Lg X = Log 10 X = Ln X / Ln 10;

 

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

 

Ln10 := 1/Ln(10);

end.

 

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

 

Lg X = Ln10 * Ln X;

 

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

 

function Lg( X : float ) : float;

Lg := Ln10*Ln(X);

 

которая позволит вычислять десятичный логарифм от 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;

SqrN := Exp( n*Ln(X) );

 

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

 

SqrtN = Exp( 1 / N * Ln(X) );

 

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

 

function SqrtN( X, n : float ) : float;

SqrtN := Exp( Ln(X) / n );

 

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

 

function Sign( X : float ): float;

if X = 0 then Sign := 1

   else Sign := X / abs(X);

 

 

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

 

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

 

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.

  

Итак, жаждите получить арксинус? Тогда их есть у меня ;o)

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

 

function ArcSin( X : float ) : float;

var y, s : float;

y := 1 – X * X;

y := abs( y );

s := sign( X );

if y = 0 then ArcSin := s * PI / 2

   else ArcSin := ArcTan( X / Sqrt( y ) );

 

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

 

function ArcCos( X : float ) : float;

var y, s : float;

y := 1 – X * X;

y := abs( y );

s := sign( X );

if X = 0 then ArcCos := PI / 2

   else ArcCos := ArcTan ( Sqrt( y ) / X );

 

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

 

function ArcCtg( X : float ) : float;

ArcCtg := PI / 2 – ArcTan( X );

 

 

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

  

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

  

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

 

При подготовке к интегрированию вычисляется коэффициент точности 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.

 

function Integral( A, B, Accuracy : float; Fx : TIntegralFunc ) : float;

const MinAccur = 0.00000000000000000001;

var T, BA, F, Fold, C, p1, p2 : float;

case SizeOf(Float) of

 6: MinAccur:=0.00000000001; {Для Real}

 4: MinAccur:=0.0000001; {Для Single}

 8: MinAccur:=0.000000000000001; {Для Double}

10: MinAccur:=0.0000000000000000001; {Для Extended}

end;

Accuracy := abs(Accuracy);

if Accuracy < MinAccur then Accuracy := MinAccur;

T := Sqr(15*Accuracy);

BA := B – A;

p1 := Fx( B ) + Fx( A );

 BA := BA / 2;

 p2 := BA + A;

 C := 0;

    C := Fx( p2 )*2 + C;

    p2 := BA*2 + p2;

 until p2 – B >= 0;

 p1 := p1 + C;

 Fold := F;

 F := ( (p1 + C) / 3 )*BA;

until Sqr( Fold – F ) –T < 0;

Integral := F;

 

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

 

Uses Math;

 

function ifunc( x : float ) : float; far;

{ где f(x[0..1]) = 1.3987174742 }

ifunc := Sqrt( 2*x + 1 );

 

function ifunc2( x : float ) : float; far;

{ где f(x[1..5]) = 0.9074539 }

ifunc2 := (x*x*x) / (x*x*x*x + 16);

 

writeln('Calc ',Integral(0,1,0,ifunc):1:10);

writeln('Calc2 ',Integral(1,5,0.01,ifunc2):1:10);

 

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

  

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

значение параметра Accuracy, что позволит получить весьма ощутимый прирост в производительности при вычислении интеграла.

  

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

 

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

 

 

 

 

 

 

Журнал > Программирование > Паскаль для новичков (Turbo Pascal, Assembler) > Паскаль для новичков (часть 24): Расширение математических возможностей Паскаля