Maxima – rozwiązywanie problemów matematycznych, część czwarta

« Część 3.: Równania. Wykresy. Granice. Pochodne. Całki. Sumy i iloczyny.

Czwarta część miała być poświęcona wykresom w zaawansowanym stopniu, ale zamiast tego będzie ona poruszała kwestię rozwiązywania równań różniczkowych i ich układów. Niniejszy tutorial jest tłumaczeniem z podręcznika „The Maxima Book”

Spis rzeczy:

Wstęp.

Uwaga.

Opis ten będzie poświęcony wyłączne równaniom różniczkowym zwyczajnym, czyli tzw. ODE.

Definiowanie ODE.

W Maximie mamy do czynienia z trzema rodzajami definiowania zwyczajnych równań różniczkowych.

Pierwszym i zarazem najprostszym sposobem jest zapis równania typu: y″ - 5y′ = -5x^2 + 2x za pomocą zwykłej pochodnej, np.: 'diff(y, x) lub 'diff(y, x, 2).

Powyższe równanie zatem byłoby wpisane do Maximy jako:

< 'diff(y, x, 2) - 5 * 'diff(y, x) = -5 * x^2 + 2 * x;
> \left( \frac{\mbox{d}^2}{\mbox{d} x^2} y \right) -5 \left( \frac{\mbox{d}}{\mbox{d} x} y \right) = -5 x^2 + 2 x

Należy zwrócić uwagę na to, że 'diff(y, x) jest poprzedzone apostrofem, co zapobiega oszacowaniu tego wyrażenia (do 0).

Drugim sposobem na definiowanie ODE, który jednocześnie pozwoli wyeliminować apostrofy, jest użycie polecenia depends, za pomocą którego informujemy Maximę, że pisząc y mamy tak na prawdę na myśli funkcję y(x). Zatem deklaracja wcześniej ustalonego przez nas równania będzie wyglądała jak następuje:

< depends
> [y(x)]
< diff(y, x, 2) - 5 * diff(y, x) = -5 * x^2 + 2 * x;
> \left( \frac{\mbox{d}^2}{\mbox{d} x^2} y \right) -5 \left( \frac{\mbox{d}}{\mbox{d} x} y \right) = -5 x^2 + 2 x

Trzecim i ostatnim sposobem jest deklaracja bezpośrednia, która polega po prostu na wpisywaniu funkcji y wprost jako zależnej od x, czyli jako y(x). Mamy więc:

< diff(y(x), x, 2) - 5 * diff(y(x), x) = -5 * x^2 + 2 * x;
> \left( \frac{\mbox{d}^2}{\mbox{d} x^2} y \right) -5 \left( \frac{\mbox{d}}{\mbox{d} x} y \right) = -5 x^2 + 2 x

Sposoby, których użyjemy zależą od poleceń jakich użyjemy do pracy ze zwyczajnymi równaniami różniczkowymi. I tak na przykład używając polecenia ode2 najczęściej będziemy stosować pierwszą i drugą formę zapisu, a do polecenia desolve będziemy zmuszeniu używać trzeciej formy.

Rozwiązywanie równań różniczkowych.

Używając ode2

Maxima potrafi rozwiązać równania różniczkowe pierwszego i drugiego rodzaju za pomocą polecenia ode2. Prototyp tego polecenia jest postaci ode2(rown, zz, zn): rozwiązywane jest równanie dane wzorem rown traktując zz jako zmienną zależną (np. y), a zn jako zmienną niezależną (np. x). Jeżeli zamiast równania rown zostanie podane wyrażenie wyr, to zostanie ono sprowadzone do postaci równania wyr = 0.

< ode2(diff(y, x, 2) - 5 * diff(y, x) = -5 * x^2 + 2 * x, y, x);
> y = k_1 \cdot e^{5 x} + \frac{x^3}{3} + k_2

W przypadku, gdy ode2 nie potrafi rozwiązać równania zwracana jest wartość FALSE.

Warunki początkowe i graniczne.

Po tym jak równanie różniczkowe zostanie rozwiązane przez ode2, możemy ustalić warunki początkowe i graniczne dla rozwiązania. Komendy pozwalające ustalić te warunki wymagają by równanie różniczkowe nie było sprecyzowane jako funkcja zmiennej, tj. diff(y, x) musi zostać użyte zamiast diff(y(x), x) by oznaczyć różniczkę.

Dla równań różniczkowych pierwszego rzędu, warunek początkowy może być ustalony za pomocą komendy ic1. Jeżeli ode2 zwróci rozwiązanie ogólne soln równania różniczkowego pierwszego rzędu, komenda ic1(soln, indvar=a, depvar=b) zwróci rozwiązanie szczególne równe b, gdy wartość zmiennej równa jest a.

< soln1:ode2(x^2*diff(y,x) + 3*x*y = sin(x)/x, y, x);
> y = \frac{\%C - \cos x}{x^3}
< ic1(soln1, x=1, y=1);
> y = - \frac{\cos x - \cos 1 - 1}{x^3}

Dla równań różniczkowych drugiego rzędu, warunki mogą być podane jako początkowe, używając komendy ic2, lub jako graniczne, używając komendy bc2. Jeżeli ode2 zwróci rozwiązanie ogólne soln równania różniczkowego drugiego rzędu, to komenda ic2(soln, indvar=a, depvar=b, diff(depvar, indvar)=c) zwróci rozwiązanie szczególne równe b i którego różniczka równa jest c, gdy wartość zmiennej jest równa a.

< eqn2: diff(y,x,2) + y = 4*x;
> \frac{\mbox{d}^2}{\mbox{d} x^2} y + y = 4x
< soln2: ode2(eqn2, y, x);
> y = \%K1 \sin x + \%K2 \cos x + 4x
< ic2(soln2, x=0, y=1, diff(y,x)=3);
> y = - \sin x + \cos x + 4x

Podobnie, jeśli ode2 zwróci rozwiązanie ogólne soln równania różniczkowego drugiego rzędu, to komenda bc2(soln, indvar=a, depvar=b, indvar=c, depvar=d) zwróci rozwiązanie szczególne równe b, gdy wartość zmiennej jest równa a i równe d, gdy wartość zmiennej jest równa c.

< bc2(soln2, x=0, y=3, x=2, y=1);
> y = - \frac{(3 \cos 2 + 7) \sin x}{\sin 2} + 3 \cos x + 4x

Metody ode2

W celu rozwiązania danego równania różniczkowego, ode2 spróbuje użyć standardowo mu dostępnych metod rozwiązywania. Metody te zostaną omówione poniżej.

Pierwszą rzeczą, którą ode2 zrobi z równaniem różniczkowym jest ustalenie czy równanie jest rzędu pierwszego czy drugiego.

Równania rzędu pierwszego.

Dla równań pierwszego rzędu, ode2 sprawdzi czy równanie pasuje do jednej z następujących kategorii i podejmie odpowiednie kroki, by równanie zostało poprawnie rozwiązane:

Liniowe równania różniczkowe niejednorodne.

Równanie różniczkowe pierwszego rzędu jest liniowe wtedy, gdy może być zapisane w postaci y′ + p(x)y = q(x). W tym przypadku, rozwiązanie jest dane jako y = (I(x) + c) / μ(x)), gdzie μ(x) jest równe eP(x) dla pierwotnej P(x) funkcji p(x), I(x) jest pierwotną μ(x)q(x), a c jest stałą.

< linode: 'diff(y,x) + x*y = x^2;
> \frac{\mbox{d}}{\mbox{d} x} y + xy = x^2
< ode2(linode,y,x);
> y = e^{-\frac{x^2}{2}}\left(\frac{i \sqrt{2}\sqrt{\pi} erf(\frac{i x}{\sqrt{2}})}{2} + xe^{\frac{x^2}{2}} + \%C \right)
Równania o zmiennych rozdzielonych.

Równanie różniczkowe pierwszego rzędu jest równaniem o zmiennych rozdzielonych, gdy może być przedstawiony w postaci M(x) = N(y) y′. W tym przypadku, rozwiązanie w postaci uwikłanej jest uzyskiwane przez całkowanie obu stron równania M(x)dx = N(y)dy.

< separableode: (3*x^2+4*x+2)=(2*y-1)*diff(y,x);
> 3x^2 + 4x + 2 = (2y - 1) (\frac{\mbox{d}}{\mbox{d} x} y)
< ode2(separableode, y, x);
> y^2 - y = x^3 + 2x^2 + 2x + \%C
Równania zupełne.

Równanie różniczkowe pierwszego stopnia jest zupełne, gdy może być przedstawione w postaci p(x,y) y′ + q(x,y) = 0, gdzie p(x,y) = ∂M(x,y)/∂y, a q(x,y) = ∂M(x,y)/∂x dla jakiegoś M(x,y). W tym przypadku rozwiązanie będzie dane w postaci uwikłanej przez M(x,y) = 0

< exactode: x^2*cos(x*y)*diff(y,x) + (sin(x*y) + x*y*cos(x*y)) = 0;
> \sin{xy} + x^2 (\frac{\mbox{d}}{\mbox{d} x} y) \cos{xy} + xy \cos{xy} = 0
< ode2(exactode,y,x);
> x \sin{xy} = \%C

Jeżeli dane równianie daje się przedstawić w postaci p(x,y)y′ + q(x,y) = 0, ale nie jest zupełne, ode2 sprawdza czy istnieje czynnik całkujący m(x,y), który sprawi, że równianie m(x,y)p(x,y)y′ + q(x,y) = 0 będzie zupełne, a w tym wypadku równanie będzie rozwiązane tak jak wyżej.

< intfactorode: (2*x*y - exp(-2*y))*diff(y,x) + y = 0;
> (2xy - e^{-2y})(\frac{\mbox{d}}{\mbox{d} x} y) + y = 0
< ode2(intfactorodee,y,x);
> x e^{2y} - \log y = \%C
Równania jednorodne.

Równianie różniczkowe pierwszego stopnia jest jednorodne, gdy może być zapisane w postaci y′ = F(y/x). W takim przypadku, podstawienie z = y/x przekształci to równanie w równanie o zmiennych rozdzielonych xz′ + z = F(z), które zostanie rozwiazane jak wcześniej.

< homode: diff(y,x) = (y/x)^2 + 2*(y/x);
> \frac{\mbox{d}}{\mbox{d} x} y = \frac{y^2}{x^2} + \frac{2y}{x}
< ode2(homode,y,x);
> -\frac{xy + x^2}{y} = \%C
Równiania Bernoullego.

Równianie postaci y′ + p(x)y = q(x)yn, n≠0,1 jest zwane równaniem Bernoullego o indeksie n. Wprowadzenie zmiennej pomocniczej zależnej z = y1-n przekształci równanie Bernoullego w równanie liniowe z′ + (1-n)p(x)z = (1-n)q(x), które zostanie rozwiązane jak wyżej.

< berode: diff(y,x) + (2/x)*y = (1/x^2)*y^3;
> \frac{\mbox{d}}{\mbox{d} x} y + \frac{2y}{x} = \frac{y^3}{x^2}
< ode2(berode,y,x);
> y = \frac{1}{\sqrt{\frac{2}{5x^5} + \%Cx^2}}
Równania rzędu drugiego

Jeżeli równanie jest drugiego rzędu, ode2 ustali czy równanie jest liniowe czy też nie. Jeżeli rówania jest liniowe, tzn. gdy moze być zapisane w postaci y″ + p(x)y′ + q(x)y = r(x), ode2 spróbuje rozwiązać równanie rozwiązując najpierw jednorodną część równaniaa, tj. y″ + p(x)y′ + q(x)y = 0. Rozwiązanie ogólne równania liniowego jednorodnego będzie postaci k1y1 + k2y2 dla stałych k1 i k2. Jeżeli r(x) ≠ 0, ode2 użyje różnych kombinacji parametrów w celu obliczenia rozwiązania szczególnego ys oryginalnego równania. Rozwiązanie ogólne równania liniowego niejednorodnego będzie więc postaci y = k1y1 + k2y2 + ys. By rozwiązać część jednorodną równania, ode2 sprawdzi czy równanie pasuje do jednej z następujących kategorii i podejmie odpowiednie kroki, by równanie zostało poprawnie rozwiązane:

Równania o stanych współczynnikach.

Jeżeli równanie ma stałe współczynniki, czyli jest postaci y″ + ay′ + by = 0, to rozwiązaniem będzie y = k1er1x + k2er2x, gdzie r1 i r2, są rozwiązaniami równania kwadratowego r2 + ar + b = 0. W przypadku gdy równanie r2 + ar + b = 0 ma podwójny pierwiastek, rozwiązaniem równania będzie y = k1erx + k2erx. W niektórych przypadkach, gdy równanie nie ma stałych współczynników, ode2 spróbuje prostej transformacji, by zredukować je do równania o stałych współczynnikach.

< ccode1: diff(y,x,2) - 3*diff(y,x) + 2*y = 0;
> \frac{\mbox{d}^2}{\mbox{d} x^2} y - 3 (\frac{\mbox{d}}{\mbox{d} x} y) + 2y = 0
< ccode2: diff(y,x,2) - 4*diff(y,x) + 4*y = 0;
> \frac{\mbox{d}^2}{\mbox{d} x^2} y - 4 (\frac{\mbox{d}}{\mbox{d} x} y) + 4y = 0
< ode2(ccode1,y,x);
> y = \%K1 e^{2x} + \%K2 e^{x}
< ode2(ccode2,y,x);
> y = (\%K2 x + \%K1) e^{2x}
Równania zupełne

Równanie różniczkowe drugiego rzędu jest zupełne, gdy może być zapisane w postaci [f(x)y′]′ + [g(x)y]′ = 0. Scałkowanie tego równania zredukuje je do równania pierwszego rzędu, które będzie rozwiązane jak wyżej.

< exactode2: x^2*diff(y,x,2) + x*diff(y,x) - y = 0;
> x^2 (\frac{\mbox{d}^2}{\mbox{d} x^2} y) + x (\frac{\mbox{d}}{\mbox{d} x} y) - y = 0
< ode2(exactode2,y,x);
> y = \%K2 x - \frac{\%K1}{2x}
Równania Eulera

Równanie postaci x2y″ + axy′ + by = 0 jest równaniem Eulera. Rozwiązanie jest dane przez y = k1xr1 + k2xr2, gdzie r1 i r2 są rozwiązaniami równania r(r-1) + ar + b = 0. W przypadku, gdy r(r-1) + ar + b = 0 ma podwójny pierwiastek, rozwiązanie jest dane przez y = k1xr + k2ln(x)xr.

< eulerode1: x^2 * diff(y,x,2) + 4*x*diff(y,x) + 2*y = 0;
> x^2 (\frac{\mbox{d}^2}{\mbox{d} x^2} y) + 4x (\frac{\mbox{d}}{\mbox{d} x} y) + 2y = 0
< elerode2: x^2 * diff(y,x,2) + 5*x*diff(y,x) + 4*y = 0;
> x^2 (\frac{\mbox{d}^2}{\mbox{d} x^2} y) + 5x (\frac{\mbox{d}}{\mbox{d} x} y) + 4y = 0
< ode2(eulerode1,y,x);
> y = \frac{\%K1}{x} + \frac{\%K2}{x^2}
< ode2(eulerode2,y,x);
> y = \frac{\%K2 \log x + \%K1}{x^2}
Uzmiennianie stałych

Jeżeli k1y1 + k2y2 jest rozwiązaniem ogólnym równania y″ + p(x)y′ + q(x)y = 0, to rozwiązanie szczególne y″ + p(x)y′ + q(x)y = r(x) może być wyznaczone gdy zmienimy stałe k1 i k2 odpowiednimi funkcjami u1 i u2 i poszukamy rozwiazania y = u1y1 + u2y2. Rozwiązanie to będzie dane, gdy u1 i u2 będą funkcjami pierwotnymi odpowiednio -y2r/(y1y2′ - y2y1′) oraz y1r/(y1y2′ - y2y1′).

< varparode: 'diff(y,x,2) + 2 * 'diff(y,x) + y = exp(x);
> \frac{\mbox{d}^2}{\mbox{d} x^2}y + 2 (\frac{\mbox{d}}{\mbox{d} x}y) + y = e^x
< ode2(varparode,y,x);
> y = \frac{e^x}{4} + (\%K2 x + \%K1) e^{-x}
Brak zmiennej zależnej

Gdy zmienna y nie występuje w równaniu różniczkowym drugiego stopnia, podstawienie ν = y′, ν′ = y″ zredukuje to równanie do równania stopnia pierwszego i może być rozwiązane jak wyżej. Gdy znajdziemy już ν, y może być wyznaczone przez scałkowanie ν.

< noyode: x * 'diff(y,x,2) + ('diff(y,x))^2 = 0;
> x (\frac{\mbox{d}^2}{\mbox{d} x^2}y) + \%DERIVATIVE^2((y,x,1)) = 0
< ode2(noyode,y,x);
> y = \int \frac{1}{\log x + \%K1}\mbox{d}x + \%K2
Brak zmiennej niezależnej

Gdy zmienna x nie występuje w równaniu różniczkowym drugiego stopnia, to robiąc zmienną zależną y tymczasowo zmienną niezależną i używając podstawienia ν = y′, ν′ = y″, zredukujemy równanie do równania stopnia pierwszego. Ponieważ różniczki są względem x, równanie to będzie mieć trzy zmienne. Może to zostać rozwiązane, ponieważ ν′ = dν/dx = (dν/dy)(dy/dx) = (dν/dy)ν oraz ν′ = dν/dx może być zastąpione przez νdν/dy. Sprawi to, że będziemy mieć do czynienia z równaniem pierwszego stopnia ze zmiennymi: zależną ν i niezależną y, które umiemy już rozwiązać używając powyższych metod. Gdy znajdziemy już ν (względem y), y będzie rozwiązaniem równania y′ = ν(y), które także umiemy już rozwiązać.

< noxode: y * 'diff(y,x,2) + ('diff(y,x))^2 = 0;
> y (\frac{\mbox{d}^2}{\mbox{d} x^2}y) + \%DERIVATIVE^2((y,x,1)) = 0
< ode2(noxode,y,x);
> \frac{y^2}{2 \%K1} = x + \%K2
Informacje o metodzie

ode2 przechowuje informacje na temat metody użytej do rozwiązania równania różniczkowego w różnych zmiennych. Zmienna method przechowuje metodę wykorzystaną do rozwiązania równania.

< ode2(separableode,y,x);
> y^2 - y = x^3 + 2 x^2 + 2 x + \%C
< method;
> SEPARABLE

W przypadku, w którym został użyty pewien czynnik, sprowadzający równanie do równania zupełnego, zmienna intfactor będzie go przechowywać.

< ode2(intfactorode,y,x);
> x e^{2y} - \log y = \%C
< method;
> EXACT
< intfactor;
> \frac{e^{2y}}{y}

Gdy rozwiązywaliśmy równanie Bernoullego lub ogólne równanie jednorodne, zmienna oindex przechowuje indeks równania.

< ode2(berode,y,x);
> y = \frac{1}{\sqrt{\frac{2}{5 x^5} + \%C} x^2}
< method;
> BERNOULLI
< odeindex;
> 3

Z kolei gdy rozwiązywaliśmy równanie niejednorodne drugiego stopnia, zmienna yp będzie przechowywać rozwiązanie szczególne pochodzące z uzmiennienia stałych.

< ode2(varparode,y,x);
> y = \frac{e^x}{4} + (\%K2 x + \%K1) e^{-x}
< method;
> VARIATIONOFPARAMETERS
< yp;
> \frac{e^x}{4}

Używając desolve

Pojedyncze równanie jednej nieznanej funkcji może być również rozwiązane za pomocą desolve. W tym przypadku nie jest konieczne podawanie ich w liście:

< de3: 'diff(f(x),x,2) + f(x) = 2 * x;
> \frac{\mbox{d}^2}{\mbox{d} x^2} f(x) + f(x) = 2 x
< desolve(de3,f(x));
> f(x) = \sin{x} \left( \frac{\mbox{d}}{\mbox{d} x} f(x) \vline_{x=0} - 2 \right) + f(0)\cos{x} + 2x

Rozwiązywanie układów równań różniczkowych

Używając desolve

Polecenie desolve

Maxima może rozwiązać układ równań różniczkowych ze stałymi współczynnikami za pomocą komendy desolve. Równania, które chcemy rozwiazać, muszą być podane za pomocą trzeciej formy (tj. diff(y(x), x)). Polecenie desolve(lista_rownan, lista_funkcji), rozwiąże układ równań danych w pierwszym parametrze (lista_rownan), gdzie lista_funkcji jest listą funkcji dla których te równania mają być rozwiązane.

< de1: diff(f(x),x) = diff(g(x),x) + sin(x);
> \frac{\mbox{d}}{\mbox{d}x} f(x) = \frac{\mbox{d}}{\mbox{d}x} g(x) + \sin{x}
< de2: diff(g(x),x) = diff(f(x),x) - cos(x);
> \frac{\mbox{d}}{\mbox{d}x} g(x) = \frac{\mbox{d}}{\mbox{d}x} f(x) - \cos{x}
< desolve([de1,de2],[f(x),g(x)]);
> \left[f(x) = e^x \left(\frac{\mbox{d}}{\mbox{d} x}g(x)\vline_{x=0} \right) - \frac{\mbox{d}}{\mbox{d} x}g(x)\vline_{x=0} + f(0) \mbox{,} g(x) = e^x \left(\frac{\mbox{d}}{\mbox{d} x}g(x)\vline_{x=0} \right) - \frac{\mbox{d}}{\mbox{d} x}g(x)\vline_{x=0} + \cos{x} + g(0) - 1\right]

Warunki początkowe

Warunki początkowe mogą być podane poleceniu desolve za pomocą polecenia atvalue. Warunki te mogą być jednak podane tylko dla 0 i muszą być podane zanim zaczniemy rozwiązywać równanie.

< atvalue(f(x),x=0,1);
> 1
< atvalue(g(x),x=0,2);
> 2
< atvalue(diff(g(x),x),x=0,3);
> 3
< desolve([de1,de2],[f(x),g(x)]);
> \left[ f(x) = 3 e^x - 2 \mbox{,} g(x) = \cos{x} + 3 e^x - 2 \right]

Metoda desolve

Algorytm desolve do rozwiązywania układów równań różniczkowych korzysta z transformaty Laplace’a. Jeśli f(t) jest zdefiniowane dla wszystkich t≥0, to wtedy tranformata Laplace’a z f jest dana przez [late bg=fff s=1x]F(s) = \mathcal{L}\left\{ f(t) \right\} = \int \limits_{0}^{\infty} e^{-st} f(t) \mbox{d}t[/latex]. Transformata Laplace’a ma ciekawą właściwość: jeżeli \mathcal{L}\left\{ f(t) \right\} = F(s), to \mathcal{L}\left\{ f'(t) \right\} = sF(s) - f(0). Transformata Laplace’a może więc zamienić układ równań różniczkowych w układ równań zwyczajnych. Jeżeli układ ten może zostać rozwiązany, to rozwiązania układu równań różniczkowych będą dane po odwróceniu transformaty Laplace’a.

Zapowiedź części piątej

W części piątej zajmiemy się wykresami (tak, wiem, część czwarta miała być temu poświęcona). Ze względów czasowych następna część pojawi się w dalekiej przyszłości, bądź wcale.

Leave a Reply