Programowanie w języku Pascal

Math, czyli matematyka studenta

  • 6 komentarzy
  • 4295 odsłon
  • Oceń ten tekst jako pierwszy
Mały przegląd funkcji TP służących do operacji na liczbach rzeczywistych,
zespolonych i macierzach - oraz kilka przydatnych przekształceń.
Wszystkie opisane procedury i funkcje znajdują się w bibliotece math
dostępnej wraz z kodem źródłowym tutaj.

Oto co nam standardowo udostępnia TP 7.0:
sin - zwraca sinus kąta podanego w radianach (2*pi = 360°);
cos - zwraca cosinus podanego kąta;
ln - logarytm naturalny (przy podstawie e) z podanej liczby;
exp - podnosi e do zadanej potęgi;
arctan - zwraca wartość kąta (w radianach!), którego tangens podajemy.

Dodatkowo jeszcze:
sqrt - pierwiastek kwadratowy (x^0.5);
sqr - kwadrat (x^2);
round - zaokrągla;
trunc - ucina część ułamkową;
int - zwraca część całkowitą (zaokrągla w dół);
frac - zwraca część po przecinku;
random - zwraca losową liczbę z przedziału <0,1>.

Niedużo tego, ale zaraz zobaczymy, co z tego można wyczarować.

Na początek zadeklarujemy sobie typ float:

type float = single;


będziemy go używać zamiast typu real, który jest wyjątkowo
nieekonomiczny jeśli chodzi o czas procesora.

Typ single jest mało dokładny (10^-7), ale za to liczenie na nim jest
sensownie szybkie. Jeśli ktoś będzie potrzebował dużej dokładności, to
zamiast single niech podstawi double albo nawet extended. Ale NIE real.

Na pierwszy ogień funkcje sin i cos:


tangens:

function tg(const x : float) : float;
var y : float;
begin
  y := cos(x);
  if y <  then tg := sin(x)/y else tg := .;
end;


Co do dzielenia przez zero to zakładam, że daje w wyniku nieskończoność.

W naszych warunkach będzie to największa dostępna na danym typie wartość. Dla single jest to 3.4e38.

Teraz cotangens:

function ctg(const x : float) : float;
var y : float;
begin
  y := sin(x);
  if y <  then tg := cos(x)/y else ctg := .;
end;


(przy okazji - jeśli mamy funkcję tg zdefiniowaną w ten sposób,
że przy dzieleniu przez zero zwraca nieskończoność, to funkcję ctg można
utworzyć tak:

function ctg(const x : float) : float;
var y : float;
begin
  y := tg(x);
  if y <  then ctg := /y else ctg := .;
end;


co jest gorszym, ale ciekawym rozwiązaniem.)

To teoretycznie wszystko w tym dziale funkcji. Kolej na arcusy:

ArcSin:
function ArcSin(const x : float) : float;
begin
  if x <  then ArcSin := ArcTan(x/sqrt (-sqr (x))) else ArcSin := .;
end;


ArcCos:
function ArcCos(const x : float) : float;
begin
  if abs(x) <  then ArcSin := ArcTan(x/sqrt (-sqr (x))) else ArcCos := .;
end;


A skoro już jesteśmy przy sinusach i im podobnych, to zdefiniujemy sobie
zestaw funkcji hiperbolicznych:

sinh - sinus hiperboliczny:

function Sinh(const x : float) : float;
begin
  sinh := (exp(x) - exp(-x))/;
end;


cosh - cosinus hiperboliczny (tzw. funkcja łańcuchowa):

function cosh(const x : float) : float;
begin
  cosh := (exp(x) + exp(-x))/;
end;


th - tangens hiperboliczny:

function th(const x : float) : float;
begin
  th := sinh(x)/cosh(x);
end;


cth - cotangens hiperboliczny:

function cth(const x : float) : float;
begin
  cth := cosh(x)/sinh(x);
end;


Jako ciekawostkę dodam, że chociaż powyższe funkcje nie mają nic
wspólnego z funkcjami sinusoidlanymi, to obowiązują dla nich
podobne prawa (np.: sinh' = cosh; cosh' = sinh; cosh2 - sinh2 = 1).

Pora na funkcje wykładnicze. Mamy funkcję exp, czyli e^x, gdzie x jest
rzeczywiste (e to liczba Eulera równa około 2.71). Czyli jak widać można
podnieść sobie e do potęgi niecałkowitej, np.: do 1/2 (czyli spierwiastkować).

A czy można podnieść dowolną liczbę rzeczywistą do takiej niecałkowitej potęgi?

No jasne! Trzeba tylko uzależnić ją od e, tzn zapisać w postaci x = e^ln(x).

Wtedy xy = (eln(x))y = e(ln(x)*y):


power - potęga rzeczywista (z liczby dodatniej!!!):

function power(const x,y : float) : float;
begin
  power := exp(ln(abs(x))*y);
end;


Jak widać x znajduje się pod logarytmem, czyli musi być dodatnie.
Natomiast y może być dowolny. Teraz można sobie sprawdzić czy to
działa: policzmy sqr(2)=4=power(2,2) i sqrt(2)=1.41=power(2,0.5).



Teraz logarytm. Turbo Pascal udostępnia nam funkcję ln, czyli
logarytm naturalny (logarytm przy podstawie e). Jeśli chcemy
policzyć algorytm przy innej podstawie, to sięgamy do wzoru z początku
liceum: loga(b) = ln(a)/ln(b).


log - logarytm przy podstawie a (a > 0 i a < 1) z liczby b (b > 0):

function log(const a,b : float) : float;
begin
  if (a > ) and (b > ) then
    log := ln(a)/ln(b)
  else
    log := .;
end;


Jeśli piszesz program, który bardzo intensywnie korzysta z którejś z
wymienionych funkcji (najczęściej jest to sin, rzadziej cos), to warto
zbudować tablicę, do niej wrzucić wartości funkcji, a potem korzystać
z tablicy pobierając przybliżone wartości funkcji. Przyspiesza to
bardzo działanie programu.

Teraz mocniejsze czary: implementacja liczb zespolonych. Nie będę dokładnie tłumaczyć, co to są liczby
zespolone, bo jeśli tego nie wiesz, to znaczy że nie i potrzebujesz tej wiedzy.

Ogólnie liczba zespolona to suma dwóch liczb rzeczywistych, z tym, że druga jest pomnożona przez tzw. jednostkę urojoną - i. Wygląda to tak: x + i*y.

Jednostka urojona ma takie miłe właściwości, że i^2 = -1 (stąd wniosek, że istnieją pierwiastki liczb ujemnych itp).

Implementacja samego typu zespolonego wygląda tak:

type
  complex = record
  case byte of
    : (re, im : float);
    : (z,  fi : float);
  end;


gdzie jednocześnie są użyte dwa typu zapisu liczb zespolonych: klasyczna
x + i*y oraz eulerowska |z|*exp(i*fi) (zawsze nieujemny moduł z, oraz
faza fi-przedział od 0 do 2pi). Tu przydadzą się funkcje do przeliczania
jednej postaci na drugą oraz do pobierania części rzeczywistej
(Realis) i urojonej (Imaginarius).

function Re(const z : complex) : float; 
begin Re := z.Re end;

function Im(const z : complex) : float; 
begin Im := z.Im end;

function zabs(const z : complex) : float; 
begin
  zabs := sqrt(sqr(z.re)+sqr(z.im));
end;

procedure efi(const z : complex;var y : complex);

begin
  y.z := zabs(z);
  if (z.re = ) and (z.im < ) then y.fi := pi/ else y.fi := abs(arctan(z.im/z.re));

  if (z.re < ) and (z.im > ) then y.fi := y.fi + pi/ else

  if (z.re < ) and (z.im <= ) then y.fi := y.fi + pi else

  if (z.re > ) and (z.im < ) then y.fi := y.fi + *pi/;

end;

procedure compl(const z : complex;var y : complex);

begin
  y.re := z.z*cos(z.fi);
  y.im := z.z*sin(z.fi);
end;


Przydałoby się też kilka procedurek obsługujących najprostsze działania na typie zespolonym:


procedure zsqr(const z : complex;var y : complex); 
begin
  y.re := sqr(z.re) - sqr(z.im);
  y.im := *z.im*z.re;
end;

procedure zmul(const a,b : complex;var y : complex); 
begin
  y.re := a.re*b.re - a.im*b.im;
  y.im := a.re*b.im + a.im*b.re;
end;

procedure zadd(const a,b : complex;var y : complex); 
begin
  y.re := a.re + b.re;
  y.im := a.im + b.im;
end;

procedure zsub(const a,b : complex;var y : complex); 
begin
  y.re := a.re - b.re;
  y.im := a.im - b.im;
end;

procedure zdiv(const a,b : complex;var y : complex); 
var
  c : complex;
begin
  c.re := b.re;
  c.im := -b.im;
  zmul(a,c,y);
  c.z := sqr(b.im) + sqr(b.re);
  y.re := y.re/c.z;
  y.im := y.im/c.z;
end;


Potrzebne też będzie pierwiastkowanie; trzeba pamiętać, że pierwiastek
stopnia n w dziedzinie zespolonej ma n rozwiązań, symetrycznie
rozłożonych wokół punktu (0,0).

procedure zsqrt(const z : complex;var x1,x2 : complex);

var
  a : complex;
begin
  efi(z,a);

  a.re := sqrt(a.re);

  x1.re := a.re*sin(a.fi);
  x1.im := a.re*cos(a.fi);

  x2.re := -x1.re;
  x2.im := -x1.im;
end;

procedure zsqrn(const z : complex;n : byte;var X : array of complex);

var
  a : complex;
  i : byte;
begin
  efi(z,a);
  a.z := power(a.z,/n);

  for i :=  to n- do
  begin
    x[i].re := a.z*cos(a.fi+ (i/n)**pi);
    x[i].im := a.z*sin(a.fi+ (i/n)**pi);
  end;
end;


Na koniec zostało jeszcze potęgowanie rzeczywiste - można je wykorzystać do
policzenia pierwszego pierwiastka zespolonego dowolnego stopnia:

procedure zpower(const z : complex;const n : float;var y : complex);
var
  a : complex;
begin
  efi(z,y);
  a.z  := power(y.z,n);
  a.fi := y.fi*n;
  y.re := a.z*cos(a.fi);
  y.im := a.z*sin(a.fi);
end;


To chyba wszystko, co najważniejsze i najniezbędniejsze w temacie liczb zespolonych.

Kolej na wektory i macierze. Najpierw definicja typu macierzowego:

const
  matrixsize =  div SizeOf(float);

type
  PVector = ^TVector;
  TVector = array[..matrixsize] of float;
  PMatrix = ^TMatrix;
  TMatrix = array[..matrixsize] of TVector; 


To dość prostackie rozwiązanie, bo nie pozwala na użycie tylko niezbędnej
ilości pamięci, więc budowa dużych macierzy jest niemożliwa (zajmują one
bardzo dużo pamięci, zbyt dużo). Częściowo można to obejść stosując wskaźniki

PVector i PMatrix, ale to połowiczne załatwienie problemu. Jednak w tym artykule
nie chodzi o oszczędność i prędkość, ale w ogóle o umożliwienie wykonywania danej operacji.

To taka mała dygresja ~:-) Teraz konkrety - podstawowe operacje na macierzach:

procedure MAdd(const A,B : TMatrix;const x,y : byte;var N : TMatrix);

var
  i,j : byte;
begin
  for i :=  to x do
  for j :=  to y do
  N[i,j] := A[i,j] + B[i,j];
end;

procedure MSub(const A,B : TMatrix;const x,y : byte;var N : TMatrix);

var
  i,j : byte;
begin
  for i :=  to x do
  for j :=  to y do
  N[i,j] := A[i,j] - B[i,j];
end;

procedure MMul(const A,B : TMatrix;const Ax,By : byte;var N : TMatrix);

var
  i,x,y : byte;

begin
  if Ax < By then exit;

  for x:= to Ax do
  for y:= to By do
  begin
    N[x,y] := ;
    for i :=  to Ax do
    N[x,y] := N[x,y] + A[i,y]*B[x,i];
  end;
end;

procedure MMulS(const A : TMatrix;const skalar : float;const x,y : byte;var N : TMatrix);

var
  i,j : byte;

begin
  for i:= to x do
  for j:= to y do
  N[j,i] := A[i,j]*skalar;
end;

procedure MTrans(const A : TMatrix;const x,y : byte;var N : TMatrix);


var
  i,j : byte;

begin
  for i:= to x do
  for j:= to y do
  N[j,i] := A[i,j];
end;

procedure MZeros(var A : TMatrix;const x,y : byte);

var
  i : byte;

begin
  for i:= to y do
  Fill16(A[,i],x*sizeof(float),);
end;


To tyle jeśli chodzi o najważniejsze rzeczy. Potęgowanie naturalne można zrealizować
poprzez pętlę mnożenia, gorzej z dzieleniem. Za jakiś czas pojawi się procedura
do odwracania macierzy (A-1), na razie jednak trzeba się jakoś obejść
bez tego. Oto przykład: rozwiązywanie układu równań A*x = b metodą Gaussa
(analogiczna do x = Ab).

function MResolve(var A : TMatrix;d : word;var x,b : TVector) : boolean;

var
  c        : float;
  n,m,i,j  : word;
  ok       : boolean;

procedure reduce(var a,n : TVector;var ab,bb : float;i : word); 
var
  y : float;
  w : word;

begin
  if a[i] =  then begin ok := false; exit; end;
  y := n[i]/a[i];
  for w := i to d do n[w] := n[w] - y*a[w];
  bb := bb - y*ab;
end;

begin
  ok := true;
  for i :=  to d- do
  for j := i+ to d do
  begin
    reduce(a[i],a[j],b[i],b[j],i);
    if not ok then break;
  end;

  if ok then
  for i := d downto  do
  begin
    if a[i,i] =  then begin ok := false; break end;
    x[i] := b[i]/a[i,i];

    for j :=  to i- do
    begin
      c := a[j,i]/a[i,i];
      b[j] := b[j] - c*b[i];
      a[j,i] := ;
    end;
  end;

  MResolve := ok;
end;


Funkcja zwróci wartość true, jeśli układ ma rozwiązanie bądź wiele rozwiązań,
jeśli układ nie ma rozwiązania, to zwrócona będzie wartość false. Po wywołaniu
funkcji w zmiennej x znajduje się wektor rozwiązań, macierz A oraz wektor b są niszczone.

Kod źródłowy biblioteki math znajduje się tutaj.

6 komentarzy

nowy121105 2006-07-23 16:33

Arcusy powinny zostac zaimplementowane w ten sposob

function ArcSin(x:float):float;
begin
if abs(x)=1 then ArcSin:=2*arctan(x) else ArcSin:=ArcTan(x/sqrt(1-sqr(x)));
end;

function ArcCos(x:float):float;
begin
ArcCos:=pi/2-ArcSin(x);
end;  

{Dla abs(x)>1  Arcusy (sinus i cosinus) (cyklometryczne) sa zespolone }

function ArcSinh(x:float):float;
begin
ArcSinh:=ln(x+sqrt(sqr(x)+1));
end;

function ArcCosh(x:float):float;
begin
ArcCosh:=ln(x+sqrt(sqr(x)-1));
end;

function ArcTanh(x:float):float;
begin
ArcTanh:=0.5*ln((1+x)/(1-x));
end;

{Dla abs(x)>1 ArcTanh jest zespolony Dla abs(x)=1 ArcTanh nie istnieje}

nowy121105 2006-07-22 18:27

Czy nie lepiej zamiast deklarowania typu complex uzyc typu string
i dopisac funkcje sprawdzajaca poprawnosc  wprowadzonych danych ?
Mozna wtedy zamiast procedur uzywac funkcji co byloby wygodniejsze.
Ciekaw jestem jak zaimplemenowac  funkcje
function re(z:string):float;
function im(z:string):float;
function complex(x,y:float):string;
uzywajac procedur str i val;
oraz
function invalidnumber(z:string):boolean
Oczywiscie nie jest to tak wygodne jak we Free Pascalu
poniewaz w TP nie mozna przeciążać operatorów ale wygodniejsze niz uzywanie
procedur

Nie dotyczy to wersji Free Pascal gdzie funkcja moze zwracac rekordy i
mozliwe jest przeciazanie a byc moze i definiowanie operatorów.

emil11 2006-06-06 16:08

super artykuł. Napewno mi się przyda. Thx

emil11 2006-06-06 16:07

super artykuł. Napewno mi się przyda. Thx

Sosna81 2005-03-21 09:33

Znalazlem blad przy mnozeniu dwóch macierzy, do poprawienia sa dwie linijki. Jest:

for i := 1 to Ax do
N[x,y] := N[x,y] + A[i,y]*B[x,i];

W takiej formie nie jest poprawny wynik, po poprawieniu na:

for i := 1 to Ay do
N[x,y] := N[x,y] + A[w,i]*B[i,y];

Jest wszystko OK.

kendas 2003-02-25 11:55

Bardzo fajny artykulik !!!