Eigenmath

(v. 137)

Portable. Freeware.

Podręcznik użytkownika

Źródło:

Eigenmath.pdf

G. Weight

Maj 2011

Przekład:

Robert Wiśniewski

http://chomikuj.pl/bobwis

Eigenmath to wieloplatformowy system algebry komputerowej napisany przez George Weigta.

Obecnie dostępny dla systemów Windows Linux i Mac OS X.

Eigenmath jest napisany w języku C, a jego kod źródłowy jest wolno dostępny.

Ten niewielki programik pozwala na wykonywanie nie tylko obliczeń numerycznych ale również symbolicznych (w  tym różniczkowanie i całkowanie), przeprowadzanie operacji macierzowych i tworzenie prostych wykresów a nawet pisanie skryptów.

W istocie rzeczy jest to zminiaturyzowana, przenośna wersja MatLab.

SPIS TREŚCI

1. Wprowadzenie

1.1. Wykładniki ujemne

1.2. Definiowanie symboli

1.3. Wykresy

1.4. Skrypty

1.5. Liczby zespolone

1.6. Algebra liniowa

2. Obliczenia

2.1. Pochodne

2.2. Gradient

2.3. Szablony funkcji

2.4. Całki

2.5. Długość łuku

2.6. Lrzywe całkowe

2.7. Pole powierzchni

2.8. Całki powierzchniowe

2.9. Teoria Greena

2.10. Teoria Stokesa

3. Dodatkowe przykłady

3.1. Wzory François Viète

3.2. Wirowość (rotacja) w postaci tensorowej

3.3. Kwantowy oscylator harmoniczny

3.4. Wodorowe funkcje falowe

3.5. Akceleratory kosmiczne Shuttle i Corvette

3.6. Stała Avogadro

3.7. Zero do potęgi zerowej

4. Funkcje wbudowane

5. Składnia

1. Wprowadzenie

Poniższe przykłady zaczerpnięte zostały z autobiograficznych notatek Vladimira Nabokova.

„Prosty samouczek objaśnił mi zbyt wcześnie istotę logarytmów i wyczytałem w (brytyjskiej publikacji dla młodzieży) o pewnym matematyku hinduskim, który dokładnie w ciągu dwóch sekund potrafił obliczać 17-ty pierwiastek np. z liczby 3529471145760275132301897342055866171392 (nie miałem pewności czy wynik jest dobry, ale sprawdziłem i wynosił on 212)”.

Możemy to sprawdzić wpisując w programie Eigenmath wyrażenie:

212^17

Po wciśnięciu klawisza ENTER, zostanie wyświetlony poniższy wynik:

3529471145760275132301897342055866171392

Tak więc Nabokow miał rację.

Możemy to dalej sprawdzić klikając przycisk Float, co prowadzi do prostszego wyniku

3.52947 x 10³⁹

Teraz możemy skorzystać z Eigenmath do obliczenia 17-tego pierwiastka z tej liczby.

N=212^17

N

N = 3529471145760275132301897342055866171392

N^(1/17)

212

Podkreśla się, że gdy do symbolu jest przypisywana w liczba, nie ukazuje się żaden wynik.

Aby zobaczyć tą wartość, trzeba wprowadzić ten symbol w wierszu wpisywania danych:

N

N = 3529471145760275132301897342055866171392

1.1. Wykładniki ujemne

Eigenmath wymaga stosowania nawiasów wokół wykładników ujemnych, np.

10^(-3)

zamiast

10^-3

Powodem tego jest to, że łączenie znaków ujemnych nie zawsze jest jednoznaczne.

Przykładowo w wyrażeniu:

X^-1/2

nie jest jasne czy wykładnikiem jest -1 czy -1/2.

Tak więc Eigenmath wymaga wprowadzania nawiasów x^(-1/2), co jest jednoznaczne.

Teraz powstaje nowe pytanie. Pomijając znak minus, czy może występować osobny znak potęgowania ^ (bez nawiasów) ?

Odpowiedź jest taka, że istotny jest tylko pierwszy znak za symbolem potęgowania. Przykładowo, poniższe wyrażenie jest traktowane jako (x¹)/2:

x^1/2

Tak więc na ogół, nawiasy sa niezbędne w wyrażeniu zawierającym wykładnik

X^(1/2)

x^1/2

1.2. Definiowanie symboli

Jak wcześniej wspomniano, Eigenmath korzysta z takiej samej składni jak stary język Fortran.

N=212^17

Przy definiowaniu symbolu, żaden wynik nie jest wyświetlany. Aby zobaczyć wynik, trzeba wprowadzić wymagany symbol:

N

N = 3529471145760275132301897342055866171392

Wszystko po pierwszej literze symbolu jest traktowane jako indeks dolny:

NA = 6.02214*10^23

N_A=6.02214*10²³

Symbol może być literą grecką:

xi=1/2

xi

Z uwagi na to, że xi jest odczytywane jako ξ, jak można wprowadzić symbol x_i ? Spróbujemy to rozwiązać w przyszłości, ale teraz xi jest odczytywane zawsze jako ξ.

Litery greckie mogą również pojawiać się w indeksach.

Amu=2.0

Amu

Aμ = 2.0

Ogólna zasada jest taka, że Eigenmath skanuje cały symbol szukając greckich liter:

alphamunu

α_μν

Spróbujemy teraz sprawdzić co się pojawi przy obliczaniu wartości symbolu. Najważniejsze jest to, że Eigenmath dokładnie oblicza podwyrażenia symboliczne.

A=B

B=C

C=D

sin(A)

sin(D)

W powyższym przykładzie, obliczanie wartości sin(A) prowadzi do wyniku sin(D), ponieważ program Eigenmath rozwiązuje symbol A zastępując go symbolem D.

Jednak wewnętrznie A jest nadal równe B. co można sprawdzić wpisując funkcję nawiasów binding:

binding(A)

B

Spróbujmy teraz wrócić na chwile do definicji symboli. Pamiętamy, że prawa strona definiowanego wyrażenia jest obliczana przed jej zamknięciem w nawiasach (przed wykonaniem funkcji binding), np.:

B=1

A=B

binding(A)

1

Funkcja binding A jest tutaj równa 1, a nie B, ponieważ B już zostało wcześniej zdefiniowane zanim pojawiało się A=B. Do poznania dokładnej wartości A stosujemy funkcję quote:

A=quote(B)

binding(A)

B

Wszystko to oznacza, że symbole mają podwójny charakter. Symbol ma wartość nawiasową binding Tora może być inna od jej obliczanej wartości, Zwykle różnica ta nie jest istotna. Funkcje quote oraz binding są tu wymienione głównie w celu zilustrowania tego co może pojawiać się wewnątrz programu. W warunkach normalnych nie trzeba korzystać z tych funkcji. Jednak jedynym istotnym wyjątkiem jest korzystanie z funkcji quote do oczyszczania symbolu.

x=3

x

x = 3

x=quote(x)

x

x

1.3. Wykresy

Funkcja draw(f,x) rysuje wykres funkcji f względem zmiennej x. Drugi argument można pominąć gdy zmienną zależną jest x lub t. Wektory xrange oraz yrange sterują skalą wykresu, np. draw(x^2).

xrange=(-1,1)

yrange=(0,2)

draw(x^2)

Funkcja Clear (lub kliknięcie przycisku Clear) resetuje zakres xrange oraz yrange przywracając go do wartości domyślnych.

Wykresy parametryczne pojawiają się gdy funkcja zwraca wektor. Funkcja trange steruje zakresem wektora.

Zakres domyślny wynosi (-π, π)

f=(cos(t),sin(t))

draw(5*f)

trange=(0,pi/2)

draw(5*f)

Poniżej przedstawiono zestaw ciekawych krzywych oraz kody do ich wykreślania.

Pierwszą z nich jest lemniskata.

clear

X=cos(t)/(1+sin(t)^2)

Y=sin(t)*cos(t)/(1+sin(t)^2)

draw(5*(X,Y))

Następną jest kardioida

r=(1-cos(t))/2

u=(cos(t),sin(t))

xrange=(-1,1)

yrange=(-1,1)

draw(r*u)

1.4. Skrypty

Poniżej przytoczono przykład rysowania wykresu linii prostej y = mx + b.

y=m*x+b

m=1/2

b=-3

draw(y)

Teraz załóżmy, że chcemy rysować ten wykres dla różnych wartości m. Możemy ponownie wszystko wpisywać, ale łatwiej jest wykonać to za pomocą skryptu. Potem wystarczy wrócić i szybko wpisywać różne wartości m oraz b tyle razy ile tylko zechcemy.

W celu utworzenia skryptu, kliknąć przycisk Edit Script. Następnie wprowadzić polecenia skryptu po jednym wierszu jak pokazano powyżej. po czym kliknąć przycisk Run Script aby zobaczyć wykres.

Program Eigenmatch uruchamia skrypt wiersz po wierszu, Każdy wiersz jest obliczany tak jak przy zwyczajnych poleceniach. Jest to kontynuowane do końca skryptu.

Po uruchomieniu skryptu można kliknąć przycisk Edit Script i w celu powrotu i dokonania zmiany skryptu.

Czasem można wydrukować skrypt np. w celu dodania do niego kilku komentarzy. Można to wykonać przez wprowadzenie wymaganego tekstu w cudzysłowach (w jednym wierszu). Przykładowo, poniższy skrypt:

„To jest liczba pi”

float(pi)

Wyświetli po kliknięciu przycisku Run Script:

To jest liczba pi

3,1415

Program Eigenmatch ma wbudowaną prostą możliwość debuggingu. Ustawienie polecenia trace na wartość 1 powoduje, że każdy wiersz skryptu będzie drukowany po uruchomieniu skryptu.

Zwykle to ustawienie powinno znajdować się w pierwszym wierszu skryptu

trace=1

Teraz po uruchomieniu skryptu drukowany będzie każdy jego wiersz.

1.5. Liczby zespolone

Gdy program Eigenmatch zostanie uruchomiony, definiuje in symbol i jako
. W pozostałych przypadkach litera i nie ma specjalnego znaczenia. Jest to zwykły symbol jaki można przedefiniować oraz korzystać z niego wg potrzeb.

Liczby zespolone można wprowadzać w postaci prostokątnej lub biegunowej:

a+i*b

a + ib

exp(i*pi/3)

Przekształcenie postaci prostokątnej w biegunową za pomocą funkcji rect upraszcza mieszane formy.

A=1+i

B=sqrt(2)*exp(i*pi/4)

A-B

rect

0

Podnoszenie do potęgi regularnych liczb zespolonych powoduje ich mnożenie.

(a+i*b)^2

Gdy a i b są liczbami a wykładnik potęgi jest ujemny, obliczenie jest wykonywane następująco:

Oczywiście, powoduje to usunięcie wielkości i z mianownika.

Poniżej podano kilka przykładów:

1/(2-i)

(-1+3i)/(2-i)

Wartość bezwzględna liczby zespolonej zwraca jej wielkość:

abs(3+4*i)

5

W związku z powyższym, można oczekiwać następującego wyniku:

abs(a+b*i)

Wynikiem nie jest,
ponieważ wymagałoby to założenia się że a i b są liczbami rzeczywistymi. Przykładowo, przyjmijmy, że a = 0, b = i. Wówczas:

oraz:

i dlatego dla pewnych a, b ∈ C, mamy

Funkcja mag jest alternatywą. Traktuje ona symbole a oraz b jako liczby rzeczywiste.

mag(a+b*i)

1.6. Algebra liniowa

Do mnożenia wektorów i macierzy stosowany jest iloczyn kropkowy (wektorowy) w postaci funkcji dot. Poniższy przykład pokazuje jak korzystać z tej funkcji oraz jak korzystać z funkcji inv do rozwiązywania układu równań liniowych AX = B względem X.

A=((3.8,7.2),(1.3,-0.9))

B=(16.5,-22.1)

X=dot(inv(A),B)

X

Co w tym jest dziwnego i dlaczego konieczne jest stosowanie funkcji dot ? Dlaczego nie korzystamy po prostu z wyrażenia
takiego jak iloczyn skalarny ? Powodem tego jest to, że program zwykle zmienia wewnętrznie kolejność czynników w celu optymalizacji procesu obliczania.

Przykładowo, symboliczna postać wyrażenia
zmieniana jest wewnętrznie na
. Z uwago na to, że iloczyn kropkowy dot nie jest przemienny, taka zmiana kolejności prowadzi do błędnego wyniku. Korzystanie z podanej wyżej funkcji do mnożenia rozwiązuje ten problem, ponieważ argument nie zmienia kolejności.

Należy podkreślić, że funkcja dot może mieć więcej niż dwa argumenty. Przykładowo, można korzystać z wyrażenia dot(A,B,C) do uzyskania iloczynu wektorowego trzech tensorów.

Poniższy przykład demonstruje relację

A=((a,b),(c,d))

inv(A)

adj(A)

det(A)

inv(A)-adj(A)/det(A)

Czasem obliczenia mogą być prostsze jeśli zostaną przeorganizowane do korzystania z funkcji adj zamiast inv. Główną ideą jest próba zapobiegania pojawiania się wyznacznika w dzielniku. Przykładowo załóżmy, że mamy macierze A i B i chcemy sprawdzić czy:

W zależności od złożoności det B, program może nie być zdolny do sprawdzenia, czy wynik jest zerem. Jeśli to ma miejsce, należy skorzystać z poniższej alternatywy:

Funkcja adj macierzy jest związana z ko faktorami w poniższy sposób:

A=((a,b),(c,d))

C=((0,0),(0,0))

C[1,1]=cofactor(A,1,1)

C[1,2]=cofactor(A,1,2)

C[2,1]=cofactor(A,2,1)

C[2,2]=cofactor(A,2,2)

C

adj(A)-transpose(C)

2. Obliczenia

2.1. Pochodne

Funkcja d(f; x) zwraca pochodną f względem zmiennej x. Zmienną x można pominąć w wyrażeniu zawierającym tylko zmienną x.

d(x^2)

2x

Poniższa tabela zestawia różne sposoby obliczania pochodnych wyższych stopni.

2.2. Gradient

Gradient f uzyskuje się przy korzystaniu z wektora x w pochodnej d(f; x).

r=sqrt(x^2+y^2)

d(r,(x,y))

Gradient f w pochodnej d(f; x) może być funkcją tensorową. Gradient wzrasta o 1 wraz ze stopniem.

F=(x+2y,3x+4y)

X=(x,y)

d(F,X)

2.3. Szablony funkcji

Funkcja f w pochodnej d(f) nie musi być zdefiniowana. Może być to szablon funkcji wyrażony poprzez jego nazwę i listę argumentów.

Program Eigenmath sprawdza listę argumentów z zadaniem do wykonania.

Przykładowo, wyrażenie d(f(x); x) jest obliczane, ponieważ f zależy od x. Wyrażenie d(f(x); y) oblicza wartość zero, ponieważ f zależy nie tylko od y.

d(f(x),x)

d(f(x),y)

0

d(f(x,y),y)

d(f(),t)

Ostatni przykład pokazuje, że pusta lista argumentów powoduje zawsze obliczanie wartości bez względu na drugi argument.

Szablony funkcji są przydatne do eksperymentowania z różnymi postaciami. Przykładowo, załóżmy że chcemy sprawdzić identyczność wyrażenia:

div(curlF) = 0

dla arbitralnej funkcji wektorowej F:

F=(F1(x,y,z),F2(x,y,z),F3(x,y,z))

curl(U)=(d(U[3],y)-d(U[2],z),d(U[1],z)-d(U[3],x),d(U[2],x)-d(U[1],y))

div(U)=d(U[1],x)+d(U[2],y)+d(U[3],z)

div(curl(F))

0

2.4. Całki

Całkowanie integral(f; x) zwraca wartość funkcji podcałkowej f względem zmiennej x.

Zmienną x można pominąć.

Można wykonywać całkowanie wielokrotne rozszerzając listę argumentów.

integral(x^2)

integral(x*y,x,y)

Funkcja defint(f, x, a, b,...) oblicza całkę oznaczoną f względem x w granicach od a do b.

Lista argumentów może być rozszerzona dla całek wielokrotnych.

Poniższy przykład oblicza całkę f=x² w dziedzinie półkola. Dla każdej wartości x wzdłuż odciętej, zakresy y wzdłuż rzędnej zmieniają się od 0 do

defint(x^2,y,0,sqrt(1-x^2),x,-1,1)

Alternatywnie można skorzystać z funkcji eval do obliczania całki oznaczonej krok po kroku.

I=integral(x^2,y)

I=eval(I,y,sqrt(1-x^2))-eval(I,y,0)

I=integral(I,x)

eval(I,x,1)-eval(I,x,-1)

Tutaj zastosowano przydatną sztuczkę. Trudniejsze całki zawierające sinus i cosinus można często rozwiązać korzystając z wykładników. Uproszczenia trygonometryczne korzystające z potęg i kątów wielokrotnych można zmienić w proste wyrażenia algebraiczne w dziedzinie wykładniczej.

Przykładem może być całka oznaczona:

którą można rozwiązać w poniższy sposób:

f=sin(t)^4-2*cos(t/2)^3*sin(t)

f=circexp(f)

defint(f,t,0,2*pi)

Poniżej podane jest sprawdzenie:

g=integral(f,t)

f-d(g,t)

0

Podstawowe prawo obliczania przedstawił James Gregory współczesny Newtonowi. Prawo to jest formalnym wyrażeniem odwrotnej relacji między całkami i pochodnymi.

Poniżej przytoczono rozwiązanie Eigenmath podstawowego prawa obliczania:

f=x

xrange=(-1,1)

yrange=xrange

draw(defint(f))

draw(integral(defint(f)))

Pierwszy wykres pokazuje, że f '(x) jest asymetryczna i dlatego łączne pole pod krzywą w granicach od -1 do 1 wynosi 0.

Drugi wykres pokazuje, że f(1) = f(-1) i dlatego dla
uzyskujemy:

2.5. Długość łuku

Załóżmy, że g(t) jest funkcją, która wykreśla krzywą. Długość łuku od g(a) do g(b) opisuje całka:

W której |g'(t)| jest długością wektora stycznego do funkcji g(t). Suma całkowa wszystkich długości łuków osiąga łączną długość od a do b. Przykładowo, spróbujemy zmierzyć długość poniższego łuku:

xrange=(0,1)

yrange=(0,1)

draw(x^2)

Odpowiednia wartość g(t) dla tego łuku wynosi:

Rozwiązanie Eigenmath ma postać:

x=t

y=t^2

g=(x,y)

defint(abs(d(g,t)),t,0,1)

float

1.47894

Zgodnie z oczekiwaniem, wynik jest większy od długości przekątnej równej
.

Wynik ten wygląda raczej na skomplikowany gdy zaczniemy od prostej paraboli.

Sprawdźmy wartości |g'(t)| aby upewnić się dlaczego tak się dzieje:

g

d(g,t)

abs(d(g,t))

Poniższy skrypt wykonuje dyskretne obliczenie długości łuku przez dzielenie krzywej na 100 części.

g(t)=(t,t^2)

h(t)=abs(g(t)-g(t-0.01))

L=0

for(k,1,100,L=L+h(k/100.0))

L

L = 1.47894

Obliczyć długość krzywej y = x^3/2 od początku układu do x = ^4/3.

x=t

y=x^(3/2)

g=(x,y)

defint(abs(d(g,x)),x,0,4/3)

Chociaż t jest podstawione w miejsce x, poprzednie rozwiązanie w istocie rzeczy nie różni się od poniższego

g=(t,t^(3/2))

defint(abs(d(g,t)),t,0,4/3)

2.6. Krzywe całkowe

Istnieją dwa rodzaje krzywych całkowych: jedne dla pól skalarnych i drugie dla pól wektorowych.

Poniższa tabela pokazuje jak są one obliczane w oparciu o długości łuków.

	Postać abstrakcyjna	Postać obliczeniowa
Długość luku
Krzywa całkowa, pole skalarne
Krzywa całkowa, pole wektorowe

Dla postaci pola wektorowego, do oznaczania jednostkowego wektora stycznego stosowany jest symbol u :

Długość wektora stycznego upraszcza się względem ds w poniższy sposób:

Po obliczeniu uzyskujemy:

oraz

gdzie C jest linią prostą od (0, 0) to (1, 1).

Jaka różnica powstaje względem pomiaru ? Pierwsza pochodna znajduje się na polem skalarnym, a druga pochodna znajduje się na polem wektorowym. Stanie się to lepiej zrozumiałe gdy uzyskamy wynik obliczeniowy:

Dlatego dla
, otrzymujemy:

x=t

y=t

g=(x,y)

defint(x*abs(d(g,t)),t,0,1)

Natomiast dla
, otrzymujemy:

x=t

y=t

g=(x,y)

F=(x,0)

defint(dot(F,d(g,t)),t,0,1)

Poniższe zadania krzywych całkowych pochodzą z podręcznika Advanced Calculus, Fifth Edition by Wilfred Kaplan.

Obliczyć całkę
wzdłuż krzywej od (0,0) do (2,2).

x=2t

y=2t

g=(x,y)

F=(y^2,0)

defint(dot(F,d(g,t)),t,0,1)

Obliczyć całkę
wzdłuż krzywej:

x=2t+1

y=t^2

z=1+t^3

g=(x,y,z)

F=(z,x,y)

defint(dot(F,d(g,t)),t,0,1)

2.7. Pole powierzchni

Niech S będzie polem powierzchni sparametryzowanej zmiennymi x i y.

Oznacza to, że
, gdzie z = f(x, y). Linie styczne w tym punkcie pola S tworzą wąski równoległobok. Powierzchnia a tego równoległoboku jest opisana wielkością iloczynu wektorowego.

Przez sumowanie wszystkich tych równoległoboków uzyskujemy łączne pole powierzchni,. A zatem:

Poniższy przykład oblicza pole powierzchni tarczy jednostkowej równoległej do płaszczyzny xy.

z=2

S=(x,y,z)

a=abs(cross(d(S,x),d(S,y)))

defint(a,y,-sqrt(1-x^2),sqrt(1-x^2),x,-1,1)

Wynikiem jest π, zgodnie z oczekiwaniem - pole powierzchni koła jednostkowego. Poniższy przykład oblicza powierzchnię z = x² + 2y jednostkowego kwadratu:

z=x^2+2y

S=(x,y,z)

a=abs(cross(d(S,x),d(S,y)))

defint(a,x,0,1,y,0,1)

Z powodów praktycznych, funkcja f(x,y) musi być bardzo prosta aby Eigenmath mógł rozwiązać całkę podwójną. Teraz obliczmy spiralę zdefiniowaną przez :

W przykładzie tym, współrzędne x, y i z są funkcjami niezależnego parametru w przestrzeni.

x=u*cos(v)

y=u*sin(v)

z=v

S=(x,y,z)

a=abs(cross(d(S,u),d(S,v)))

defint(a,u,0,1,v,0,3pi)

float

10.8177

2.8. Całki powierzchniowe

Całka powierzchniowa jest podobna do łapania wiatru w żagle.

Inaczej mówiąc, obliczamy poniższą całkę:

w której F⋅n jest ilością wiatru padającego na wąski, prostopadły równoległobok dA. Całka ta sumuje po całej powierzchni żagla.

Niech S będzie powierzchnią żagla sparametryzowaną przez x i y (w modelu tym kierunek z jest zgodny z wiatrem).

Na podstawie właściwości iloczynu wektorowego uzyskujemy poniższe wyrażenie dla jednostki normalnej n oraz dla dA:

Wynika stąd, że:

Przykładowo, obliczmy poniższa całkę powierzchniową:

gdzie:
, S jest powierzchnią,
, a n jest górą .

Podkreśla się, że powierzchnie te przecinają się w płaszczyźnie koła xy.

Korzystając z reguły prawej dłoni przecinamy x oraz y, co prowadzi do górnego punktu n.

Poniższy kod Eigenmatch oblicza tą całkę powierzchniową.

Symbole f i h stosowane są tutaj jako zmienne tymczasowe.

z=1-x^2-y^2

F=(x*y^2*z,-2*x^3,y*z^2)

S=(x,y,z)

f=dot(F,cross(d(S,x),d(S,y)))

h=sqrt(1-x^2)

defint(f,y,-h,h,x,-1,1)

2.9. Teoria Greena

Teoria Greena mówi nam, że:

Przykład 1: Obliczyć
wokół koła
korzystając z teorii Greena. Przy próbie obliczania stwierdzamy, że Eigenmatch nie może rozwiązać bezpośrednio całki podwójnej względem x i y.

Dlatego korzystamy zamiast tego ze współrzędnych biegunowych.

P=2x^3-y^3

Q=x^3+y^3

f=d(Q,x)-d(P,y)

x=r*cos(theta)

y=r*sin(theta)

defint(f*r,r,0,1,theta,0,2pi)

Funkcja podcałkowa defint ma postać f*r, ponieważ
.

Teraz spróbujemy obliczyć boczną krzywą całkową wg teorii Greena i sprawdźmy czy uzyskamy ten sam wynik.

Musimy tu zastosować pewną sztuczkę przekształcając sinus i cosinus na wykładniki, aby Eigenmath mógł znaleźć wymagane rozwiązanie.

x=cos(t)

y=sin(t)

P=2x^3-y^3

Q=x^3+y^3

f=P*d(x,t)+Q*d(y,t)

f=circexp(f)

defint(f,t,0,2pi)

Przykład 1: Obliczyć wg teorii Greena obie strony
względem tarczy
. Najpierw obliczamy krzywą całkową wzdłuż granic tarczy. Podkreśla się, że promień tarczy wynosi 2.

--Line integral

P=1-y

Q=x

x=2*cos(t)

y=2*sin(t)

defint(P*d(x,t)+Q*d(y,t),t,0,2pi)

--Surface integral

x=quote(x) --remove parametrization of x

y=quote(y) --remove parametrization of y

h=sqrt(4-x^2)

defint(d(Q,x)-d(P,y),y,-h,h,x,-2,2)

--Bonus point: Compute the surface integral using polar coordinates.

f=d(Q,x)-d(P,y) --do before change of coordinates

x=r*cos(theta)

y=r*sin(theta)

defint(f*r,r,0,2,theta,0,2pi)

defint(f*r,theta,0,2pi,r,0,2) --try integrating over theta first

W tym przypadku Eigenmath rozwiązuje obydwie formy całki biegunowej. Jednak w przypadkach, w których Eigenmath nie może uzyskać rozwiązania całki podwójnej, można zmienić kolejność całkowania.

2.10. Teoria Stokesa

Teoria Stokesa oferuje poniższe równoważniki krzywych całkowych i całek powierzchniowych.

gdzie F = (P, Q, R).

Dla wielkości S sparametryzowanej względem x i y, uzyskujemy:

W wielu przypadkach, przekształcenie całki zgodnie z teorią Stokesa może rozwiązać problem sprowadzając go do łatwiejszej postaci.

Załóżmy, że F = (x, y, z) oraz że S jest częścią paraboloidy
, która znajduje się nad płaszczyzną xy.

Obwód tej paraboloidy jest kołem
.

Obliczyć krzywe całkowe oraz całki powierzchniowe.

Tu również stwierdzimy, że treba skorzystać w tym celu ze współrzędnych biegunowych aby uzyskać sukces.

--Surface integral

z=4-x^2-y^2

F=(y,z,x)

S=(x,y,z)

f=dot(curl(F),cross(d(S,x),d(S,y)))

x=r*cos(theta)

y=r*sin(theta)

defint(f*r,r,0,2,theta,0,2pi)

--Line integral

x=2*cos(t)

y=2*sin(t)

z=4-x^2-y^2

P=y

Q=z

R=x

f=P*d(x,t)+Q*d(y,t)+R*d(z,t)

f=circexp(f)

defint(f,t,0,2pi)

3. Dodatkowe przykłady

3.1.Wzory François Viète

François Viète był pierwszym, który odkrył dokładną formułę obliczania wartości liczby π.

Poniżej pokazano tą formułę:

Załóżmy, że
oraz
. Możemy więc napisać:

Rozwiązując to wyrażenie względem π, uzyskujemy:

Teraz skorzystajmy z Eigenmatch do obliczenia π zgodnie formułą Viète. Oczywiście, nie możemy obliczać wszystkich elementów aż do nieskończoności. Zamiast tego ograniczymy się do pierwszych 9-ciu składników, co zapewnia nam uzyskanie dokładności do 6 cyfr znaczących.

a(n)=test(n=0,0,sqrt(2+a(n-1)))

float(2*product(k,1,9,2/a(k)))

Funkcja a(n) może być wywoływana n-razy, a więc łącznie istnieje tu 54 wywołań a(n). Korzystając z innego algorytmu opartego na szablonach zmiennych możemy uzyskać ten wynik w 9 krokach.

a=0

b=2

for(k,1,9,a=sqrt(2+a),b=b*2/a)

3.2. Wirowość (rotacja) w postaci tensorowej

Funkcja wirowości Curl wektora może być wyrażona w postaci wektorowej w poniższy sposób:

gdzie
jest tensorem Levi-Civita.

Poniższy skrypt pokazuje, ze formuła ta jest równoważna obliczaniu wirowości starym sposobem:

-- Define epsilon.

epsilon=zero(3,3,3)

epsilon[1,2,3]=1

epsilon[2,3,1]=1

epsilon[3,1,2]=1

epsilon[3,2,1]=-1

epsilon[1,3,2]=-1

epsilon[2,1,3]=-1

-- F is a generic vector function.

F=(FX(),FY(),FZ())

-- A is the curl of F.

A=outer(epsilon,d(F,(x,y,z)))

A=contract(A,3,4) --sum across k

A=contract(A,2,3) --sum across j

-- B is the curl of F computed the old fashioned way.

BX=d(F[3],y)-d(F[2],z)

BY=d(F[1],z)-d(F[3],x)

BZ=d(F[2],x)-d(F[1],y)

B=(BX,BY,BZ)

-- Are A and B equal? Subtract to find out.

A-B

Poniżej pokazano odmianę poprzedniego skryptu. Iloczyn
jest obliczany za pomocą jednego wiersza kodu. Ponadto, drugi iloczyn oraz kontrakcja wzdłuż k jest teraz obliczana za pomocą iloczynu skalarnego.

F=(FX(),FY(),FZ())

epsilon=zero(3,3,3)

epsilon[1,2,3]=1

epsilon[2,3,1]=1

epsilon[3,1,2]=1

epsilon[3,2,1]=-1

epsilon[1,3,2]=-1

epsilon[2,1,3]=-1

A=contract(dot(epsilon,d(F,(x,y,z))),2,3)

BX=d(F[3],y)-d(F[2],z)

BY=d(F[1],z)-d(F[3],x)

BZ=d(F[2],x)-d(F[1],y)

B=(BX,BY,BZ)

-- A

3.3. Kwantowy oscylator harmoniczny

Łączną energię E, energię kinetyczną K oraz energię potencjalną V wyraża równanie:

E = K + V

Odpowiadająca formuła kwantowego oscylatora harmonicznego ma poniższą postać:

gdzie n jest liczbą całkowitą reprezentującą kwantyfikowane wartości energii.

Powyższe równanie ma następujące rozwiązanie:

w którym Hn(x) jest wielomianem Hermite'a względem x.

Poniższy kod Eigenmath sprawdza E = K + V dla n = 7.

n=7

psi=exp(-x^2/2)*hermite(x,n)

E=(2*n+1)*psi

K=-d(psi,x,x)

V=x^2*psi

E-K-V

3.4. Wodorowe funkcje falowe

Funkcje falowe wodoru
są rozwiązaniami równania różniczkowego:

gdzie n jest liczbą całkowita reprezentującą kwantyfikację łącznej energii, natomiast r jest odległością promieniową elektronu, Operator Laplasjan we współrzędnych sferycznych ma postać:

Ogólną postacią
jest:

gdzie L jest wielomianem Laguerre, P jest wielomianem Legendre, natomiast l oraz m są liczbami całkowitymi spełniającymi poniższe zależności:

Ogólną postać można wyrazić jako iloczyn promieniowej funkcji R i harmoniki sferycznej Y.

Poniższy kod Eigenmath sprawdza E = K + V dla n , l. M = 7; 3; 1.

laplacian(f)=1/r^2*d(r^2*d(f,r),r)+

1/(r^2*sin(theta))*d(sin(theta)*d(f,theta),theta)+

1/(r*sin(theta))^2*d(f,phi,phi)

n=7

l=3

m=1

R=r^l*exp(-r/n)*laguerre(2*r/n,n-l-1,2*l+1)

Y=legendre(cos(theta),l,abs(m))*exp(i*m*phi)

psi=R*Y

E=psi/n^2

K=laplacian(psi)

V=2*psi/r

simplify(E-K-V)

0

3.5. Akceleratory kosmiczne Shuttle i Corvette

Wahadłowce kosmiczne Shuttle uzyskują przyspieszenia od 0 do 17 000 mil/godzinę w ciągu 8 minut. Statki kosmiczne Corvette uzyskują przyspieszenia od 0 do 60 000 mil/godzinę w ciągu 4,5 sekund. Poniższy skrypt porównuje te dwa wyniki:

vs=17000*"mile"/"hr"

ts=8*"min"/(60*"min"/"hr")

as=vs/ts

as

vc=60*"mile"/"hr"

tc=4.5*"sec"/(3600*"sec"/"hr")

ac=vc/tc

ac

"Time for Corvette to reach orbital velocity:"

vs/ac

vs/ac*60*"min"/"hr"

Poniżej pokazano te wyniki uzyskane po uruchomieniu skryptu. Wynika stąd, że wahadłowiec uzyskuje przyspieszenie dwukrotnie większe niż statek kosmiczny.

Czas osiągnięcia prędkości orbitalnej przez statek kosmiczny wynosi:

3.6. Stała Avogadro

Zaproponowano zdefiniowanie stałej Avogadro dokładnie na poziomie 84446886 do potęgi trzeciej . Liczba ta odpowiada idealnemu sześcianowi atomów zawierającemu 84,446,886 wzdłuż każdej krawędzi. Można sprawdzić, że różnica między wartością zaproponowaną i zmierzoną wynosi zaledwie 6.0221415 ± 0.0000010) x 10²³ atomów.

A=84446886^3

B=6.0221415*10^23

A-B

0.0000010*10^23

Jak widać, zaproponowana wartość mieści się w granicach błędu eksperymentu. Możemy rozłożyć na czynniki zaproponowaną liczbę:

factor(A)

3.7. Zero do potęgi zerowej

Poniższy przykład tworzy wykres funkcji
. Wykres ten pokazuje dlaczego konwencja zapisu 0⁰ = 1 ma sens.

f(x)=abs(x^x)

xrange=(-2,2)

yrange=(-2,2)

draw(f)

Z wykresu tego widać jak wyrażenie 0⁰ = 1 tworzy linię ciągłą w okolicy x = 0. Teraz możemy zobaczyć jak zachowuje się wyrażenie x^x na płaszczyźnie zespolonej:

f(t)=(real(t^t),imag(t^t))

xrange=(-2,2)

yrange=(-2,2)

trange=(-4,2)

draw(f)

4. Funkcje wbudowane

Funkcja	Opis	Przykłady
abs(x)	Zwraca wartość bezwzględną lub długość wektora x. Dla liczb zespolonych można stosować funkcję mag(z)	P(x,y) abs(P)
adj(m)	Zwraca dopełnienie macierzy m
and(a,b,...)	Zwraca logiczne AND określanych wyrażeń
arccos(x)	Zwraca odwrotność cosinusa x
arccosh(x)	Zwraca odwrotność cosinusa hiperbolicznego x
arcsin(x)	Zwraca odwrotność sinusa x
arcsinh(x)	Zwraca odwrotność sinusa hiperbolicznego x
arctan(x)	Zwraca odwrotność tangensa x
arctanh(x)	Zwraca odwrotność tangensa hiperbolicznego x
arg(z)	Zwraca kąt liczby zespolonej z
ceiling(x)	Zwraca najmniejszą liczbę całkowitą nie mniejszą od x
check(x)	W skrypcie, gdy określana wartość x jest prawdziwa, kontynuuje obliczanie. W przeciwnym razie zatrzymuje.
choose (n,k)	Zwraca wartość wyrażenia
circexp(x)	Zwraca wyrażenie x z funkcją kołową przekształconą w postać wykładniczą. Czasem to upraszcza wyrażenie.
coeff(p,x,n)	Zwraca współczynnik x^n wielomianu p
cofactor(m,i,j)	Zwraca kofaktor macierzy m względem wiersza i oraz kolumny j.
con(j)	Zwraca sprzężenie liczby zespolonej z
contract(a,i,j)	Zwraca tensor a sumowany po wskaźnikach i oraz j. Gdy i oraz j są pominięte, wówczas stosowane są indeksy 1 oraz 2. Funkcja contract(m) jest równoważna ścieżce macierzy m.
cos(x)	Zwraca cosinus x
cosh(x)	Zwraca cosinus hiperboliczny x
cross(u,w)	Zwraca iloczyn wektorowy wektorów u oraz v
curl(u)	Zwraca wirowość (rotację) wektora u
d(f,x)	Zwraca pochodną funkcji f względem zmiennej x
defint(f,x,a,b)	Zwraca całkę oznaczoną funkcji f względem zmiennej x w granicach od a do b. Listę argumentów można rozszerzać na całki wielokrotne. Przykładowo: defint(f,x,a,b,y,c,d)
deg(p,x)	Zwraca stopień wielomianu p względem x
denominator(x)	Zwraca mianownik wyrażenia x
det(m)	Zwraca wyznacznik macierzy m
do(a,b,…)	Oblicza listę argumentów od lewej do prawej. Zwraca wynik dla ostatniego argumentu
dot(a,b,…)	Zwraca iloczyn skalarny tensorów
draw(f,x)	Wykreśla funkcję f względem zmiennej x
erf(x)	Zwraca funkcję błędu względem x
erfc(x)	Zwraca komplementarną funkcję błędu względem x
eval(f,x,n)	Zwraca wartość f obliczoną dla x = n
exp(x)	Zwraca wartość e^x
expand(r,x)	Zwraca cząstkowe rozszerzenie ułamkowe stosunku dwóch wielomianów r względem x	expand(1/(x^3+x^2),x
expcos(x)	Zwraca cosinus x w postaci wykładniczej	expcos(x)
epsin(x)	Zwraca sinus x w postaci wykładniczej	epsin(x)
factor(n)	Rozkłada liczbę całkowitą n na czynniki pierwsze	factor (12345)
factor(p,x)	Rozkłada wielomian n na czynniki względem x. Ostatni argument może być pominięty. Lista argumentów może być rozszerzona na wielomiany kilku zmiennych. Przykładowo, factor(p,x,y) rozkłada wielomian p na czynniki x, a następnie na czynniki y.	factor (125*x*3-1)
Factorial(n)	Zwraca silnię liczby całkowitej n. Alternatywnie, można stosować polecenie n!	10!
filter(f,a,b)	Zwraca wyrażenie f z usuwanymi członami a, b, itd.	1/a + 1/b + 1/c filter(last,a)
float(x)	Przekształca x na liczbę zmiennoprzecinkową	sum(n,0,20,(-1/2)^n) float(last)
floor(x)	Zwraca największa liczbę całkowitą, nie większą od x
for(i,j,k,a,b,…)	Dla i = j poprzez k oblicza a, b, itp.	x=0 y=2 for(k,1,9,x=sqrt(2+x), y=2*y/x) float(y)
gcd(a,b….)	Zwraca największy wspólny dzielnik
hermite(x,n)	Zwraca wielomian Hermite'a n-tego stopnia względem x
hilbert(n)	Zwraca macierz Hilberta rzędu-n
imag(z)	Zwraca część urojoną liczby zespolonej z
inner(a,b,...)	Zwraca iloczyn wewnętrzny tensorów. Funkcja ta spełnia taką samą rolę jak iloczyn skalarny
inegral(f,x)	Zwraca całkę f względem zmiennej x
inv(m)	Zwraca odwrotność macierzy m
isprime(n)	Zwraca wartość 1 gdy n jest liczbą pierwszą. W przeciwnym razie zwraca 0.	Isprime(2^53-111)
laguerre(x,n,a)	Zwraca wielomian Laguerre n-tego stopnia. Gdy pominiemy parametr a, stosowane jest a = 0.
lcm(a,b,...)	Zwraca najmniejszą wspólną wielokrotność
leading(p,x)	Zwraca wiodący współczynnik wielomianu p względem x
legendre(x,n,m)	Zwraca wielomian Legendre n-tego stopnia. Gdy pominiemy parametr m, stosowane jest m = 0.
log(x)	Zwraca logarytm naturalny x
mag(z)	Zwraca rozmiar liczby zespolonej z
mod(a,b)	Zwraca resztę z dzielenia a przez b
not(x)	Neguje wynik określanej wartości x
nroots(p,x)	Zwraca wszystkie pierwiastki (rzeczywiste i zespolone) wielomianu p względem x. Pierwiastki te obliczane są w sposób numeryczny. Współczynniki wielomianu p mogą być rzeczywiste lub zespolone.
numerator(x)	Zwraca licznik wyrażenia x
or(a,b,...)	Zwraca logiczne OR określanej wartości x
outer(a,b,...)	Zwraca iloczyn zewnętrzny tensorów
polar(z)	Przekształca liczbę zespoloną z do postaci biegunowej
prime(n)	Zwraca n-tą liczbę pierwszą
print(a,b,...)	Oblicza wyrażenie i drukuje wynik. Przydatne w pętli for
product(i,j,k,f)	Zwraca wartość
quote(x)	Zwraca wyrażenie x bez obliczania
quotient(p,q,x)	Zwraca iloraz wielomianów względem x
rank(a)	Zwraca liczbę indeksów w tensorze a. Skalar nie ma indeksów, a więc jego rząd jest równy 0
rationalize(x)	Wstawia wszystko do wspólnego mianownika	rationalize(a/b+b/a)
real(z)	Zwraca część rzeczywistą liczby zespolonej z
rect(z)	Zwraca liczbę zespoloną z w postaci prostokątnej
roots(p,x)	Zwraca takie wartości x, że wielomian p(x) = 0. Wielomian musi by rozkładany na czynniki całkowite
simplify(x)	Zwraca x w najprostszej (uproszczonej) postaci
sin(x)	Zwraca sinus x
sinh(x)	Zwraca sinus hiperboliczny x
sqrt(x)	Zwraca pierwiastek kwadratowy x
stop	W skrypcie zatrzymuje obliczenia
subst(a,b,c)	Podstawia a za b w wyrażeniu c i zwraca wynik
sum(i,j,k,f)	Zwraca wartość
tan(x)	Zwraca tangens x
tanh(x)	Zwraca tangens hiperboliczny x
taylor(f,x,n,a)	Zwraca rozwinięcie Taylora f względem x wokół punktu a. Argument n jest stopniem rozwinięcia. Gdy pominiemy a, stosowane jest a = 0	taylor(1/cos(x),x,4)
test(a,b,c,d)	Gdy a jest prawdą, wówczas zwracane jest b, ponadto gdy c jest prawdą, wówczas zwracane jest d, itd. Gdy liczba argumentów jest nieparzysta, wówczas zwracany jest ostatni argument, gdy wszystkie pozostałe są fałszywe
transpose(a,i,j)	Zwraca transpozycję tensora względem jego wskaźników i oraz j. Gdy pominiemy i oraz j, wówczas stosowane są wartości 1 oraz 2. Dlatego można wykonać transpozycję macierzy zawierającej 1 argument	A=((a,b),(c,d)) transpose(A)
unit(n)	Zwraca macierz jednostkową n x n	unit(2)
zero(i,j,…)	Zwraca tensor zerowy o wymiarach i, j. itp. Funkcja przydatna do tworzenia tensora i ustawiania wartości jego składników

5. Składnia

Matematyka	Eigenmath	Inna postać lub komentarz
-a	-a
a + b	a+b
a - b	a-b
ab	a*b	a b ze spacją między nimi
	a/b
	a/b/c
a^2	a^2
	a^(1/2)	sqrt(a)
	a^(-1/2)	1/sqrt(a)
a(b + c)	a*(b+c)	a (b+c) ze spacją między nimi
f(a)	f(a)
	(a,b,c)
	((a,b),(c,d))
T^12	T[1.2]	Dostęp do składników tensora
2 km	2*”km”	Jednostki miary w cudzysłowach

Tekst komentarzy nie może zawierać polskich znaków diakrytycznych (przypis tłumacza).

Williamson and Trotter, Multivariable Mathematics, p. 598.

Kaplan, Advanced Calculus, p. 313

Wilfred Kaplan, Advanced Calculus, 5th Edition, 287.

Fox, Ronald and Theodore Hill. \An Exact Value for Avogadro's Number." American Scientist 95 (2007) 104-107.

Proponowana liczba w artykule ma aktualnie wartość 84446888³.

W dodatku, autor zredukował ja do 84446886³ aby była podzielna przez 12.

See www.physorg.com/news109595312.html

- 33 -

---