Politechnika Warszawska - Wydział Samochodów i Maszyn Roboczych
Teoria Maszyn i Podstawy Automatyki
plik pomocniczy do zajęć projektowych
Sebastian Korczak
10.12.2018

/opracowane w wxMaxima 17.10.1/

Parametry liczbowe definiujemy z użyciem dwukropka.
Mozemy uzywac liczba całkowitych, z rozwinięciem dziesiętnym i ułamków.
Polecenia oddzielamy średnikiem. Enterem przechodzimy do nowej linii.
Kod uruchamiamy kombinacja CTRL+ENTER (lub enterem klawiatury numerycznej).
Po uruchomieniu kodu jest on przesyłany do głównej aplikacji, wykonywane są obliczenia
a następnie zwracane wyniki opatrzone oznaczeniem zmiennej w nawiasie.
Maxima stara sie używać w pierwszej kolejności obliczeń symbolicznych i ułamków.
Jeśli chcemy otrzymać wartość z rozwinięciem dziesiętnym używamy funkcji float().

a: 4;
b: 1.5;
c: 1/3;
d: 3-2/5;
z: float(d);

Na liczbach i parametrach możemy wykonywać dowolne działania.

b+a*d;
sqrt(b)-d;
a+g-%omega^2*sqrt(z);
sqrt(2);

Możemy wykorzystać liczbę pi albo e, jak również jednostkę urojoną - oznaczamy je procentem

%pi;
%e;
%i;

ZASTOSOWANIE WXMAXIMA W AUTOMATYCE


Transmitancje operatorowe proponuję deklarować w następujący sposób:

G1: 5;
G2: 1/(2*s+3);
G3: k/(3.5*s^2+4.1*s+2);

Na transmitancjach zawierajacych zarowno liczby jak i parametry mozna wykonywać różne
operacje matematyczne.

Gz: (G1*G3)/(1+G1*G3*G2);

Widzimy, że maxima wstawia wzory i zwraca postać bez żadnych uproszczeń.
Jeśli potrzebujemy uproszczonej formy zapisu możemy użyć dwóch funkcji:
expand() - rozbicie znalezionych iloczynow nawiasow na sumy
ratsimp() - uproszczenie wyrazenia do postaci wielomianowej

Gz1: expand(Gz);
Gz2: ratsimp(Gz);

Maxima czasem informuje nas, że dokonuje jakichś liczbowych uproszczeń,
które były niezbędne w uproszczeniach, np. zastępuje ułamki poprzez liczby
zmiennoprzecinkowe (rat: replaced...).

TRANSFORMATA LAPLACE'A

Transformate Laplace'a wyrażenia liczymy funkcją
laplace(funkcja czasu, zmienna czasu, zmienna zaspolona)

x: a;
X: laplace(x,t,s);
y: %e^(-2*t);
Y: laplace(y,t,s);

Odwrotną transformatę Laplace'a wyrażenia liczymy funkcją
ilt(funkcja zespolona, zmienna zespolona, zmienna czasu)
(skrót od Inverse Laplace Transform)

Z: 2/(3*s+2);
z: ilt(Z,s,t);
U: 2/(2*s^2+3*s+2);
u: ilt(U,s,t);
F: 2*s/(s^3+4*s^2+5*s+2);
f: ilt(F,s,t);

Jesli program nie jest w stanie policzyc odwrotnej transformaty laplacea, to znaczy, ze
mianownik transmitancji posiada pierwiastki zespolone, a z takimi program sobie nie radzi.
Zwraca wtedy polecenie, którego nie mógł rozwiązać.

R: 1/(s^3+4*s^2+11*s+2);
r: ilt(R,s,t);

Pierwiastki wielomianu mozemy jednak policzyć sami analitycznie z pomoca funkcji
solve (wielomian, zmienna), może ona zawierać nawet parametry. Przykłady:

M1: s^3+4*s^2+5*s+2;
solve(M1,s);
M2: s^2+k*s+2;
solve(M2,s);
M3: s^2+2*s+8;
solve(M3,s);
W: aa*s^2+bb*s+cc;
W: solve(W,s);

Przypominam, że jednostka urojona w maximie oznaczana jest przez %i, a w automatyce używamy j

Jeśli metoda analityczna sobie nie radzi, to możemy policzyc w sposob numeryczny (przyblizony) funkcją
 allroots(wielomian, zmienna),
przy czym wielomian moze miec tylko jedną zmienną

M4: s^3+4*s^2+5*s+2;
allroots(M4,s);
M5: s^4+4*s^3+4*s^2+5*s+2;
allroots(M5,s);

Tworząc transmitancję widmową możemy uzyć funkcji podstawiania
subst(co wstawić, zamiast czego, gdzie)

Gz: 2/(2*s^2+3*s+2);
Gzw: subst(%i*w, s, Gz);

Część rzeczywistą i urojoną z transmitancji widmowej możemy otrzymać poleceniami:

P: realpart(Gzw);
Q: imagpart(Gzw);

W przypadku bardziej skomplikowanej transmitancji dodatkowo możemy wspomagać
się funkcjami upraszczającymi

Gz: (s^3+8)/(5*s^6+3*s^5+2*s^4+s^3+s^2+s+4);
Gzw: subst(%i*w, s, Gz);
P: realpart(Gzw);
Q: imagpart(Gzw);

RYSOWANIE WYKRESÓW

Wykres funkcji jednej zmiennej tworzymy poleceniem wxplot2d(....), gdzie:
[t,0,15] - oznacza zakres zmiennej na osi x
grid2d - rysuje siatkę
ylabel i xlabel - definiuja podpisy osi
legend false lub true - ukrywa lub wyswietla legedę

q: sin(3*t);
wxplot2d(q, [t, 0, 15],grid2d,[ylabel,"y(t)"],[xlabel,"t"],[legend,false])$

Jeśli zależy nam na modyfikacji wielkości wykresu można użyć analogicznej funkcji plot2d,
która otworzy wykres w drugim oknie z ikonami do edycji

plot2d(q, [t, 0, 15],grid2d,[ylabel,"y(t)"],[xlabel,"t"],[legend,false])$

Wzmocnienie i przesunięcie fazowe najszybiej liczy sie tak:
(uwaga, w maximie jest dostepna tylko funkcja logarytmu naturalnego log(), aby wyliczyc logarytm 
 dziesietny z liczby należy logarytm naturalny tej liczby podzielić przez logarytm naturalny z 10)

A: ratsimp(sqrt(P^2+Q^2));
L: (20*log(A)/log(10));
Fi: atan2(Q,P);

Wykres Bodego rysujemy z zastosowaniem logarytmicznej skali na osi poziomej [logx].
Odpowiednio dobieramy zakres częstości większych od zera.

Bode_wzm: wxplot2d(L, [w, 0.01, 100],grid2d,[ylabel,"L(w)[dB]"],[xlabel,"w[rad/s] (log)"],[legend,false],[logx])$
Bode_faza: wxplot2d(Fi, [w, 0.01, 100],grid2d,[ylabel,"Fi(w)[rad]"],[xlabel,"w[rad/s](log)"],[legend,false],[logx])$

Charakterystyka NYQUISTA wymaga narysowania wykresu parametrycznego.
Musi więc przebiegać według poniższego wzoru
(uwaga! mozna uzyc tylko zmiennej t)

XX:2/(t^2+4)
$ YY:-t/(t^2+4) $
wxplot2d([parametric, XX, YY,[t, 0.001, 1000],[nticks, 100]],[xlabel,"P(w)"],[ylabel,"Q(w)"],grid2d)$

Uwaga! rysowanie funkcji parametrycznej wymaga użycia zmiennej t (błąd programu, który nie zezwala na inną zmienną).
Możemy więc wcześniej liczone funkcje P i Q skorygować zamieniając "w" na "t" i dopiero rysując

P1: subst(t, w, P);
Q1: subst(t, w, Q);
wxplot2d([parametric, P1, Q1,[t, 0.01, 10],[nticks, 100]],
    [xlabel,"P1"],[ylabel,"Q1"],[legend,false],
    [style,[lines,2]],[gnuplot_preamble, "set grid;"], grid2d)$

Liczenie zapasu modułu:
1. szukamy częstości w1 dla której Q(w)=0
2. Liczymy wartość P dla tej omegi w1
3. Zapas modułu wyniesie 1+P(w1)

w1: find_root(Q=0, w, 0.01, 1000);
Pw1: subst(w1,w,P);
dM: 1+Pw1;
rad: float(%pi+subst(w1,w,Fi));
deg: float(%*180./%pi);

Liczenie zapasu fazy:
1. szukamy częstości w1 dla której A(w)=1
2. Liczymy wartość Fi dla tej omegi w1
3. Zapas modułu wyniesie PI+P(w1)

w1: find_root(A=1, w, 0.01, 1000);
FIw1: subst(w1,w,Fi);
dFIrad: float(%pi+FIw1);
FIw1deg: float(%*180./%pi);

Macierze definiujemy poleceniem matrix podajac po przecinku wiersze macierzy. Weznacznik liczy funkcja determinant(...)

D: matrix( [ 5, 2, 1 ],
                 [ 3, 2, 1 ],
                [ 2, 2, 2 ]  );
determinant(D);

Rozkład na ułamki proste jest możliwy z użyciem funkcji partfrac(...)
(ale działa tylko dla rzeczywistych pierwiastków mianownika)

partfrac(2*s/(2*s^2+4*s+2),s);