Maxima
Z Wikibooks, biblioteki wolnych podręczników.
Przykładowe zadania z rozwiązaniami w programie Maxima (wxMaxima wersja 0.7.1, Maxima 5.13.0, GNU Common Lisp GCL 2.6.8)
Spis treści |
[edytuj] Wprowadzanie danych
Każde wprowadzone polecenie kończymy średnikiem ; bądź symbolem dolara $
Zakończenie polecenia średnikiem powoduje, że Maxima oblicza wyrażenie i wyświetla wyniku jego działania:
(%i1) 2*2; (%o1) 4
Zakończenie polecenia symbolem dolara $ powoduje, że Maxima oblicza wyrażenie, przechowuje jego wynik (w tym przykładzie zapisany w %o2) lecz nie wyświetla go:
(%i2) 2*2$ (%i3)
Wynik jest dostępny zmienną % (ostatni obliczony wynik)
(%i3)%; (%o3) 4
bądź, w naszym przypadku, %o2
(%i4)%o2; (%o4) 4
Wyjątki : Wewnątrz instrukcji block polecenia kończy się przecinkiem, oprócz ostatniego, które nie ma swojego końcowego znaku.
[edytuj] Przebieg programu
Interpreter Maximy przyjmuje:
- pojedyncze proste polecenie (1 linia)
- polecenie złożone (wiele linii)
- programy, czyli pliki wsadowe zawierające poleceniami. Są to pliki tekstowe (ANSI) z rozszerzenie mac. Otwieramy np. przez Menu/File/Batch file.
/* zawartość pliku z rozszerzeniem mac */ i:1; block ( loop, i:i+1, if i<100 then go(loop) ); i;
[edytuj] Liczby
[edytuj] Liczby zespolone
Postać algebraiczna (kanoniczna) c:x+y*%i
(%i1) c:5+3*%i; (%o1) 3*%i+5
część rzeczywista:
(%i2) x:realpart(c); (%o2) 5
część urojona:
(%i3) y:imagpart(c); (%o3) 3
Moduł :
(%i4) r:float(abs(c)); (%o4) 5.830951894845301;
Argument :
(%i5) t:float(carg(c)); (%o5) 0.54041950027058;
Postać wykładnicza:
(%i10) c1:r*%e^(%i*t); (%o10) sqrt(34)*%e^(%i*atan(3/5))
Powrót do formy algebraicznej :
(%i11) rectform(%); (%o11) 3*%i+5
Postać trygonometryczna :
(%i32) c2:r*(cos(t)+%i*sin(t)); (%o32) 5.830951894845301*((3*%i)/sqrt(34)+5/sqrt(34))
Powrót do formy algebraicznej :
(%i33) rectform(%); (%o33) (17.4928556845359*%i)/sqrt(34)+29.1547594742265/sqrt(34) (%i34) float(%), numer; (%o34) 3.0*%i+4.999999999999999
Porównywanie dwóch liczb zespolonych:
(%i12) IsEqual(c1,c2,eps):=
if abs(realpart(c1)-realpart(c2))<=eps
and
abs(imagpart(c1)-imagpart(c2))<=eps
then true
else false;
Sprawdzamy działanie funkcji:
(%i13) IsEqual(1.1, 1.27, 0.1); (%o13) false
Głowna wartość argumentu liczby zespolonej w obrotach :
carg_t(z):= block( [t], t:carg(z)/(2*%pi), /* teraz jest w obrotach ale w zakresie (-pi,pi] */ if t<0 then t:t+1, /* zmiana zakresu z (-pi,pi] do [0, 2*pi) */ return(t) )$
Pierwiastek kwadratowy
csqrt(z):= block( [t,w,re,im], t:abs(z)+realpart(z), if t>0 then (re:sqrt(t/2), im:imagpart(z)/sqrt(2*t)) else (im:abs(z), re:0), w:re+im*%i, return(w) )$
[edytuj] Funkcje
Definicja prostej funkcji:
(%i1) f(z,c):=z*z+c; (%o1) f(z,c):=z*z+c
Jej pierwszej pochodnej względem zmiennej z:
(%i2) m(z,c):=diff(f(z,c),z,1); (%o2) m(z,c):=diff(f(z,c),z,1)
Sprawdzamy wzór pochodnej :
(%i3) m(z,c); (%o3) 2*z
Definicja funkcji rekurencyjnej:
(%i12) F(n, z, c) :=
if n=1 then f(z,c)
else f(F(n-1, z, c),c);
(%o12) F(n,z,c):=if n=1 then f(z,c) else f(F(n-1,c,z),c)
Sprawdzamy wzór:
(%i13) F(2,z,c); (%o13) (z^2+c)^2+c
[edytuj] Rozwiązywanie równań
[edytuj] Pojedyncze równanie
Chcemy znaleźć punkty okresowe dla funkcji Fc(z)=z*z+c dla c=0+1*i (z i c są liczbami zespolonymi).
W tym celu najpierw określamy wartość c:
(%i1) c:1*%i; (%o1) %i
[edytuj] Okres=1 (punkty stałe)
Rozwiązujemy:
(%i7) p:float(rectform(solve([z*z+c=z],[z]))); (%o7) [z=0.62481053384383*%i-0.30024259022012,z=1.30024259022012-0.62481053384383*%i]
Szukamy punktu odpychającego (źródło, repeller), czyli beta :
if (abs(2*rhs(p[1]))<1)
then block(
beta:rhs(p[1]),
alfa:rhs(p[2])
)
else block(
alfa:rhs(p[1]),
beta:rhs(p[2])
);
Pokazujemy punkt stały typu alfa:
(%i20) alfa; (%o20) 0.62481053384383*%i-0.30024259022012
Pokazujemy punkt stały typu beta:
(%i21) beta; (%o21) 1.30024259022012-0.62481053384383*%i
Pokazujemy, że suma alfa i beta jest równa 1
(%i10) p:solve([z*z+c=z], [z]); (%o10) [z=-(sqrt(1-4*%i)-1)/2,z=(sqrt(1-4*%i)+1)/2] (%i14) s:radcan(rhs(p[1]+p[2])); (%o14) 1
Rysujemy punkty:
(%i15) xx:makelist (realpart(rhs(p[i])), i, 1, length(p)); (%o15) [-0.30024259022012,1.30024259022012] (%i16) yy:makelist (imagpart(rhs(p[i])), i, 1, length(p)); (%o16) [0.62481053384383,-0.62481053384383] (%i18) plot2d([discrete,xx,yy], [style, points]);
Można dodać też punkt z=1/2 (alfa i beta są symetryczna względem 1/2):
(%i31) xx:cons(1/2,xx); (%o31) [1/2,-0.30024259022012,1.30024259022012] (%i34) yy:cons(0,yy); (%o34) [0,0.62481053384383,-0.62481053384383] (%i35) plot2d([discrete,xx,yy], [style, points]);
[edytuj] Punkty o okresie = 2
(%i2) solve([(z*z+c)^2+c=z], [z]); (%o2) [z=-(sqrt(-4*c-3)+1)/2,z=(sqrt(-4*c-3)-1)/2,z=-(sqrt(1-4*c)-1)/2,z=(sqrt(1-4*c)+1)/2]
Pokazujemy że z1+z2 = -1
(%i4) s:radcan(rhs(p[1]+p[2])); (%o4) -1
Pokazujemy że z2+z3=1
(%i6) s:radcan(rhs(p[3]+p[4])); (%o6) 1
[edytuj] Punkty o dowolnym okresie
(%i41) giveRoots(period,c):=allroots(%i*F(period,z,c)); (%o41) giveRoots(period,c):=allroots(%i*F(period,z,c))
[edytuj] Układ równań
[edytuj] liniowych
[edytuj] nieliniowych
Układ równań nieliniowych rozwiązujemy za pomocą komendy algsys (przykłady Stavrosa Macrakisa). Wynik zależy co uznamy za parametr a co za zmienną.
(%i2) algsys([a*x+y=3,x+b*y=1],[x,y]); (%o2) x=(3*b-1)/(a*b-1),y=(a-3)/(a*b-1)
(%i3) algsys([a*x+y=3,x+b*y=1],[a,b]); (%o3) a=-(y-3)/x,b=-(x-1)/y
(%i4) algsys([a*x+y=3,x+b*y=1],[x,b]); (%o4) x=-(y-3)/a,b=(y+a-3)/(a*y)
(%i5) algsys([a*x+y=3,x+b*y=1],[y,a]); (%o5) y=-(x-1)/b,a=(x+3*b-1)/(b*x)
Jeśli liczba równań jest mniejsza niż liczba zmiennych to wynikiem będzie równanie parametryczne
(%i1) algsys([a*x+y=3,x+b*y=1],[a,b,x,y]); (%o1) a=(3*%r1-%r1*%r2)/(%r1*(1-%r1*%r2)),b=%r1,x=1-%r1*%r2,y=%r2
[edytuj] Przykład 1 :
Rozwiązywanie układu równań: Definujemy pierwsze równanie:
(%i1) z*z+c=z; (%o1) z^2+c=z
Definiujemy drugie równanie ( t jest to kąt w obrotach ; 1 obrót = 360 stopni) :
(%i2) 2*z=%e^(2*t*%pi*%i); (%o2) 2*z=%e^(2*%i*%pi*t)
Ponieważ 2 równanie ma tylko 1 zmienną eliminujemy zmienną z pierwszego równania :
(%i3) eliminate([%i1,%i2],[z]); (%o3) [%e^(4*%i*%pi*t)-2*%e^(2*%i*%pi*t)+4*c]
Rozwiązujemy równanie %o3 względem zmiennej c :
(%i4) solve([%e^(4*%i*%pi*t)-2*%e^(2*%i*%pi*t)+4*c], [c]); (%o4) [c=-(%e^(4*%i*%pi*t)-2*%e^(2*%i*%pi*t))/4]
Definiujemy funkcję :
(%i34) define(c1(t),rhs(s[1])); (%o34) c1(t):=-(%e^(4*%i*%pi*t)-2*%e^(2*%i*%pi*t))/4
Sprawdzamy działanie :
(%i35) c1(0); (%o35) 1/4
Rysujemy wykres ( główny komponent zbioru Mandelbrota ) :
plot2d ([parametric, realpart(c1(t)), imagpart(c1(t)), [t,0,2*%pi],[nticks,2000]])$
[edytuj] Rysowanie
[edytuj] Rysowanie 1 funkcji
Definicja funkcji f1:
f1(x):=a*x*(1-x)$
Rysowanie funkcji y:=f1(x) w zakresie 0 <= x <=1:
plot2d(f1(x),[x,0,1]);
[edytuj] Rysowanie kilku funkcji na 1 wykresie
Definicja funkcji f1:
f1(x):=%pi/2;
Definicja funkcji f2:
f2(x):=-%pi/2;
Rysowanie:
plot2d([f1(x),f2(x),atan(x)],[x,-5,5],[plot_format, gnuplot]);
Inny przykład:
plot2d([round(x),x,x*x,floor(x),ceiling(x),1/x,log(x)], [x,0,5],[y,0,5],[plot_format, gnuplot])$
[edytuj] plot2D
Dla niektórych funkcji f(x) typowa instrukcja:
plot2d (f(x), [x, 0, 1])
działa prawidłowo, dla innych daje błąd[1].
Rozwiązaniem jest wywołanie funkcji w postaci : '(f(x)) w liście argumentów plot2d. Przykłady:
f(x) := if x > 1 then 2/(x^2+1) else x^2$
plot2d ('(f(x)), [x, 0, 3])$
g(x) := quad_qag (t^2, t, 0, x, 2) [1]$
plot2d ('(g(x)), [x, -2, 2])$
W tych przykładach, jeśli użyjemy f(x) and g(x) bez znaku ' to pojawi się błąd.
Proszę zwrócić uwagę :
'(f(x))
nie
'f(x).
[edytuj] Rysowanie funkcji w różnych postaciach
Rysujemy koło o promieniu =1 i środku = 0
- w postaci biegunowej:
load(draw); draw2d(polar(1,theta,0,2*%pi));
- w postaci parametrycznej:
plot2d ([parametric, cos(t), sin(t), [t,0,2*%pi],[nticks,80]])$
- w postaci uwikłanej
load(implicit_plot); z:x+y*%i; implicit_plot (abs(z) = 1, [x, -4, 4], [y, -4, 4]);
Rysujemy lemniscatę
- w postaci uwikłanej :
load(implicit_plot); z:x+y*%i; f(n):=abs(z^n-1); implicit_plot (f(4) = 1, [x, -2, 2], [y, -2, 2]);
[edytuj] Rysowanie zbioru punktów
- tworzymy 2 listy :
- xx z wartościami x
- yy z wartościami yy
- wyświetalmy listy używając procedury draw2d z pakietu draw
load(draw); draw2d( points(xx,yy));
[edytuj] Rysowanie funkcji 2 zmiennych
plot3d(atan2(y,x), [x,-5,5], [y,-5,5], [plot_format,gnuplot])$
[edytuj] Zapisywanie wykresu do pliku eps
W przypadku Windows plik będzie w katalogu c:\
plot2d([f1(x),f2(x)],[x,0,1],[gnuplot_term,ps], [gnuplot_out_file,"log.eps"]);
[edytuj] Biblioteki graficzne
[edytuj] Różnice między Maximą a innymi językami
[edytuj] Instrukcja for
W Maximie dozwolone jest użycie wartości niecałkowitych jako indeksów pętli for :
for x:0.5 step 0.001 thru 1.0 do
