Maxima

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)

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 wynik 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.

Przebieg programu

Program rysujący zbiór Julia

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,
• moduły z definicjami funkcji ( package ) napisane w Maximie (*.mac) lub Lispie ( *.lisp) za pomocą funkcji load lub z Menu ( load package )

Przykład prostego polecenia, nadajemy zmiennej a wartość 2 :

a:2;

Jak widać Maxima wykorzystuje znak dwukropka dla nadania wartości

Przykład pliku wsadowego ( możemy go zapisać jako np. a.mac ) :

/* zawartość pliku wsadowego zawierającego program z rozszerzeniem mac */ 
i:1;
block
(
  loop,
  i:i+1,
  if i<100 then go(loop)
);
i;

Przykłady ładowania modułu :

(%i2) load("complex_dynamics.lisp");
(%o2) C:/PROGRA~1/MAXIMA~1.2/share/maxima/5.19.2/share/dynamics/complex_dynamics.lisp

lub

(%i1) load("dynamics"); 
(%o1) C:/PROGRA~1/MAXIMA~1.2/share/maxima/5.19.2/share/dynamics/dynamics.mac

Funkcje możemy skompilować. Wtedy szybciej działają :

compile(all);

Liczby

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

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)
)$

Sortowanie listy liczb zespolonych wg argumentu :

l2_ordered: sort(l2, lambda([z1,z2], is(carg(z1) < carg(z2))))$

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

Zmienne lokalne funkcji

Pamiętaj że zmienne lokalne musisz zadeklarować przed użyciem ( porównaj konstrukcje w Lispie ).

give_circle(center,radius,n_points):=block(
[t,dt,x,y,xy], /* deklaracja */
t:0, 
dt:(2*%pi)/n_points,
x:realpart(center)+radius*cos(t),
y:imagpart(center)+radius*sin(t),
xy:makelist ([x,y], k, 1, 1),
for i:1 thru (n_points-1) step 1 do
( t:t+dt,
 x:realpart(center)+radius*cos(t),
 y:imagpart(center)+radius*sin(t),
 xy:cons([x,y],xy)
 ),
return(xy)
);

Tak używamy tej funkcji :

center:0;
radius:1;
n_points:10;
xy:give_circle(center,radius,n_points);
load(draw);
draw2d( points(xy))$

Porównaj ten kod z :

n_points:1000;
dt:(2*%pi)/n_points;
t:0*%pi;
/* compute coordinates or root point */ 
x:cos(t);
y:sin(t);
/* save coordinates to list*/ 
xy:makelist ([x,y], i, 1, 1);
for i:1 thru (n_points-1) step 1 do
( t:t+dt,
  x:cos(t),
  y:sin(t),
  xy:cons([x,y],xy)
);
load(draw);
draw2d( points(xy))$

Rozwiązywanie równań

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

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]);

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

Punkty o dowolnym okresie

(%i41) giveRoots(period,c):=allroots(%i*F(period,z,c));
(%o41) giveRoots(period,c):=allroots(%i*F(period,z,c))

Metoda Newtona

load (newton1);

definiujemy równanie :

eq1:z^3-1;

szukamy pojedynczego rozwiązania w okolicy %i :

s:newton(eq1,z,%i,0.1);

otrzymujemy :

 (%o19) -(0.33333333333333*((-(0.33333333333333*((%i+0.33333333333333*(-%i-1))^3-1))/(%i+0.33333333333333*(-%i-1))^2+%i
+0.33333333333333*(-%i-1))^3-1)) /(-(0.33333333333333*((%i+0.33333333333333*(-%i-1))^3-1))/(%i+0.33333333333333*(-%i-1))^2+%i+0.33333333333333*(-%i-1))^2
-(0.33333333333333*((%i+0.33333333333333*(-%i-1))^3-1))/(%i+0.33333333333333*(-%i-1))^2+%i+0.33333333333333*(-%i-1)

zmieniamy formę wyniku :

s1:rectform(float(s));

otrzymujemy :

0.86816551188735*%i-0.50879080328932

Układ równań

liniowych

$\begin{cases} 2x-3=0\\ x-y=0 \end{cases}$

$(x,y)=(\frac{3}{2},\frac{3}{2})$

solve([2x-3=0,x-y=0],[x,y])

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

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]])$

Rysowanie

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]);

inny sposób :

load(draw);
draw2d(explicit(x,x,-1,1));

Topologczny sinus :

plot2d (sin(1/x), [x, 0, 1])$

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]);

wykrs funkcji atan

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])$

plot2D

Dla niektórych funkcji f(x) typowa instrukcja:

plot2d (f(x), [x, 0, 1])

działa prawidłowo, dla innych daje błąd.

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).

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,-%pi,%pi],[nticks,80]], [x, -4/3, 4/3])$

• 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(2) = 1, [x, -2, 2], [y, -2, 2]);

• w postaci parametrycznej

plot2d ([parametric, cos(t), sin(t)*cos(t), [t,-%pi,%pi],[nticks,80]], [x, -4/3, 4/3])$

Rysowanie zbioru punktów

Rysowanie punktów 2D. Grafika i kod żródłowy

• 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));

albo tworzymy 1 listę xy z parami wartości x i y i rysujemy :

load(draw);
draw2d( points(xy));

Topologiczny grzebień . Grafika i kod żródłowy

Topologiczny grzebień :

load(draw); /* by Mario Rodríguez Riotorto  http://riotorto.users.sourceforge.net/gnuplot/index.html  */

draw2d(
 title     = "Topologist\\47s comb", /* Syntax for postscript enhanced option; character 47 = '  */
 xrange        = [0.0,1.2],
 yrange        = [0,1.2],
 file_name = "comb2",
 terminal      = svg,
 points_joined = impulses, /* vertical segments are drawn from points to the x-axis  */
 color         = "blue",
 points(makelist(1/x,x,1,1000),makelist(1,y,1,1000)) )$

Rysowanie funkcji 2 zmiennych

plot3d(atan2(y,x), [x,-5,5], [y,-5,5], [plot_format,gnuplot])$

wykres funkcji atan2

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"]);

Biblioteki graficzne

• Draw

Tworzenie pliku SVG

BeginSVG(file_name,cm_width,cm_height,i_width,i_height):=
block(
 destination : openw (file_name),
 printf(destination, "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>~%"),
 printf(destination, "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"~%\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">~%"),
 printf(destination,"<svg width=\"~d cm\" height=\"~d cm\" viewBox=\"0 0 ~d ~d\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">~%",
 cm_width,cm_height,i_width,i_height),
 return(destination)
);
CircleSVG(dest,center_x,center_y,_radius):=printf(dest,"<circle cx=\"~f\" cy=\"~f\" r=\"~f\" fill=\"white\" stroke=\"black\" stroke-width=\"2\"/>~%",
 center_x,center_y,_radius);
CloseSVG(destination):=
 (
 printf(destination,"</svg>~%"),
 close (destination)
 );
/* ---------------------------------------------------- */
cmWidth:10;
cmHeight:10;
iWidth:800;
iHeight:600;
radius:200;
centerX:400;
centerY:300;
f_name:"b.svg";
/* ------------------------------------------------------*/ 
f:BeginSVG(f_name,cmWidth,cmHeight,iWidth,iHeight);
CircleSVG(f,centerX,centerY,radius); 
CloseSVG(f);

Różnice między Maximą a innymi językami

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

Zmienne lokalne funkcji

Zmienne lokalne musisz zadeklarować przed użyciem ( porównaj konstrukcje w Lispie ).

Źródła

• | Podstawowy kurs Maximy autorstwa Rafała Topolnickiego
• | Matematyka wyższa z Maximą autorstwa Wojciecha Młocka
• | Przewodnik po Maximie
• | Macierze w programie wxMaxima