http://sulfurzona.com/
News
Service
Magazine
Software (Battle City Game, Wallpaper manager, Superpad, VG-NOW, Puzzle Game, Netler Internet Browser, ..)
Dune Game (Dune III, Dune IV, Cheats, Forum, ..)
Games free
Turbo Pascal (Assembler, Docs, Sources, Debbugers, ..)
Books (Docs for developers)
Forum
Guest book
Компьютерная диагностика двигателя автомобиля (адаптер К-линии)Компьютерная диагностика двигателя автомобиля (адаптер К-линии)
 
 
 

Паскаль для новичков (часть 24)

 

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

 

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

 
Паскаль для новичковОднажды, когда мне понадобилось реализовать сложный алгоритм с использованием логарифмов, возведений в произвольную степень, вычисление определённых интегралов, я был неприятно удивлён, что среди встроенных алгебраических функций среды 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 по основанию 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;
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.
  
Итак, жаждите получить арксинус? Тогда их есть у меня ;o)
Вызывайте функцию 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
 
При подготовке к интегрированию вычисляется коэффициент точности 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;
begin
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 );
repeat
 BA := BA / 2;
 p2 := BA + A;
 C := 0;
 repeat
    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;
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, что позволит получить весьма ощутимый прирост в производительности при вычислении интеграла.
  
Чтобы сократить время вычисления при достаточно высокой точности, следует использовать более сложные методы численного интегрирования, как методы Бодэ, Ньютона-Котеса, Уэддля, Чебышева, Гаусса, и др.
 
Пришёл, увидел, а побеждать уже нечего. ;o)
 
Литература
1. Дьяконов В.П. Справочник по расчётам на микрокалькуляторах – М.:Наука, 1989. – 462 с.
 
Продолжение следует…
 
© Владислав Демьянишин
 
 
Вы находитесь на официальном сайте Владислава Демьянишина - разработчика игры Dune IV (Dune 4). На нашем сайте Вы можете бесплатно скачать игры Dune IV (Dune 4), Battle City (Танчики с Dendy/Nintendo), читы к играм и многое другое. Также Вы можете скачать бесплатно программы и полезные утилиты. Все программы чистые, т.е. не содержат вирусов и иного вредоносного ПО.
 
Среди доступных программ есть мобильная читалка книг, менеджер переноса файлов с фото- и видеокамер на компьютер, текстовый редактор, WYSIWYG редактор, 3D аниматор, GIF аниматор, AVI аниматор, пакетный конвертор изображений, редактор электрических схем, программа для скриншотов, диспетчер тем рабочего стола и другие.
 
На нашем сайте можно не только бесплатно скачать игры, но и документацию и книги по программированию на MIDLetPascal, Turbo Pascal 6, Turbo Pascal 7, Borland Pascal, по программированию устройств Sound Blaster, Adlib, VESA BIOS, справочник Norton Guide и много другой полезной информации для программистов, включая примеры решения реальных задач по созданию резидентных программ. Предлагаю также посетить Марья искусница - сайт о рукоделии (http://mariya-iskusnica.ru).
 
 

Журнал > Программирование > Паскаль для новичков (Turbo Pascal, Assembler) > Паскаль для новичков (часть 24): Расширение математических возможностей Паскаля
 
 
 
108
 
ВКонтакте
Facebook
 
 
 
На главную страницу На предыдущую страницу На начало страницы
 
 
Украинский портАл Украина онлайн Рейтинг@Mail.ru Рейтинг Сайтов YandeG Rambler's Top100