Rozwiązywanie liniowych równań różniczkowych jednorodnych drugiego rzędu metodą Frobeniusa

Projekt zaliczeniowy kursu: Mathematica II: aplikacje
Paweł Nowak
Fizyka doświadczalna III rok
Wydział Fizyki, Astronomii i Informtyki Stosowanej
Uniwersytet Jagielloński
Kraków 2010

Zawartość projektu

1. Zawartość projektu 			s2
2. Wstęp					s3
    2.1. Trochę teorii 				s3
3. Rozwiązania					s4
4. Podsumowanie                                         	s5
5. Literatura					s6

Wstęp

Jednorodne linowe równania rożniczkowe zwyczajne o zmiennych współczynnikach często pojawiają się w różnego rodzaju zagadnieniach przyrodniczych. Bardzo rzadko można znaleźć rozwiązanie ścisłe, dlatego wypracowano wiele metod ich rozwiązywania. Jednym ze sposobów jest metoda szeregów potęgowych, określana metodą Frebeniusa. Za pomocą tego nootbooka można poszukiwać rozwiązań w postaci szeregów. Opis metody Frebeniusa; założenia i zakres stosowalności jest przedstawiony na kolejnym slajdzie.

Trochę teorii

Załóżmy, że równanie różniczkowe możemy sprowadzić do postaci:

(1.a)             	P (x) y''(x) + Q (x) y'(x) + R (x) y(x) = 0

lub poniższego kształtu:

(1.b)             	y''(x) + p (x) y'(x) + q(x) y(x) = 0

gdzie: 
		p(x) = ((Q(x))/(P(x))) i q(x) = ((Q(x))/(R(x))).
		
Rozwiązanie ogólne równania różniczkowego (1.b) ma postać:

		y(x) = (Subscript[c, 1])(Subscript[y, 1])(x) +(Subscript[c, 2]) (Subscript[y, 2])(x).
		
Tutaj (Subscript[y, 1])(x) i (Subscript[y, 2])(x) są dwoma liniowo niezależnymi szczególnymi rozwiązania równania (1.b).

Zakładamy że funkcje (x-(Subscript[x, 0]))p(x) i (((x-(Subscript[x, 0])))^2)p(x) są analityczne dla |x-(Subscript[x, 0])| < r (promień zbieżności) i są rozwijalne w szeregi

(2)                   (x-(Subscript[x, 0]))p(x) =  (∑_(n=0)^∞)(Subscript[p, n])(((x-(Subscript[x, 0])))^n),

(3)		 (((x -(Subscript[x, 0])))^2)q(x) = (∑_(n=0)^∞)(((Subscript[q, n])(x-(Subscript[x, 0])))^n).

Postulujemy rozwiązanie postaci:

(4)		y(x)=(∑_(n=0)^∞)(Subscript[a, n])(((x-(Subscript[x, 0])))^(n+s)) z (Subscript[a, 0])≠0.

Podstawiając (4) do (1.b) otrzymamy równanie wskaźnikowe:

(5)                  s(s-1)+(Subscript[p, 0])s+(Subscript[q, 0])=0.

Rozwiązując równanie (5) otrzymujemy dwa rozwiązania (Subscript[s, 1]) i (Subscript[s, 2]); należy rozpatrzyć dwa przypadki:

1) kiedy  (Subscript[s, 1]) - (Subscript[s, 2]) ∉ ℕ
2) kiedy  (Subscript[s, 1]) -(Subscript[s, 2]) ϵ ℕ

Dla przypadku 1) rozwiązaniami liniowo niezależnymi są:

 (Subscript[y, 1])(x)=(∑_(n=0)^∞)(Subscript[a, n])(((x-(Subscript[x, 0])))^(n+(Subscript[s, 1])))

i

 (Subscript[y, 2])(x)=(∑_(n=0)^∞)(Subscript[b, n])(((x-(Subscript[x, 0])))^(n+(Subscript[s, 2]))).

Natomiast dla przypadku 2) mamy dwa liniowo niezależne rozwiązania w postaci:

 (Subscript[y, 1])(x)=(∑_(n=0)^∞)(Subscript[a, n])(((x-(Subscript[x, 0])))^(n+(Subscript[s, 1])))

oraz

 (Subscript[y, 2])(x)=(∑_(n=0)^∞)(Subscript[b, n])(((x-(Subscript[x, 0])))^(n+(Subscript[s, 2]))) + c (Subscript[y, 1])(x) ln(x - (Subscript[x, 0]))

zaś stała c może być równa zeru.

Mając już pierwsze szczególne rozwiązanie (Subscript[y, 1])(x) możemy skonstruować drugie niezależne rozwiązanie wg wzoru:

(Subscript[y, 2])(x)= d (Subscript[y, 1])(x)∫((Exp(-∫p(x)dx))/((Subscript[y, 1])(((x))^2)))dx,
gdzie d jest dowolną stałą.

uwagi

Należy zauważyć, że przypadek gdy (Subscript[s, 1]) = (Subscript[s, 2 ])jest przypadkiem szczególnym 2) dlatego nie jest rozpatrywany oddzielnie.

Rozwiązania

Teraz czas na praktyczną stronę metody.

ClearAll["Global`*"]

Twoje równanie jest zapewne w postaci: operator L(y, x) = 0. 
Wykonaj poniższą komórkę i wpisz operator L(y, x) w zadanej formie.

Operator3[y_,x_]=Input["Podaj operator w formie: P(x)y''[x] + Q(x)y'[x] + R(x)y[x]",P[x]y''[x] + Q[x]y'[x] + R[x]y[x]]

Teraz określ ilość wyrazów znaczących w rozwinięciu w szereg wokół punktu (Subscript[x, 0]) = 0.

n=Input["Podaj wyraz rozwinięcia"]

P[x_]=Coefficient[Operator3[y,x],y''[x]];
Q[x_]=Coefficient[Operator3[y,x],y'[x]];
R[x_]=Coefficient[Operator3[y,x],y[x]];
Print["P[x]=",P[x]]
Print["Q[x]=",Q[x]]
Print["R[x]=",R[x]]

p[x_]=((Q[x])/(P[x]));
q[x_]=((R[x])/(P[x]));
Print["p[x]=",p[x]]
Print["q[x]=",q[x]]

Policzmy teraz współczynniki i wstawmy je do równania wskaźnikowego.

wspolczynniki = Limit [{x p[x],x^2q[x]},x->0];
rownaniewskaznikowe=r(r-1)+wspolczynniki[[1]]r+wspolczynniki[[2]]//Simplify;
wskazniki=Solve[rownaniewskaznikowe==0,r];
rr=Sort[r/.wskazniki,#1>#2&]
l1:=rr[[1]]
l2:=rr[[2]]
dl:=l1-l2

Teraz pora na pierwsze szczególne rozwiązanie

y[x_]:=x^l1(∑_(i=0)^n)a[i]x^i+O[x]^(n+1)//Simplify;
ukladrownan=LogicalExpand[Operator3[y,x]==0];
wspolczynikirownania=Table[a[i],{i,1,n}];
wartosciwspolczynnikow=Solve[ukladrownan,wspolczynikirownania];
y1[x_]=y[x]/.wartosciwspolczynnikow/.a[0]->1//Simplify//First;
yr1[x_]=Normal[y1[x]/.a[0]->1];

Sprawdzamy pierwsze rozwiązanie.

Operator3[y1,x]

Drugie szczególne rozwiązanie

yp1[x_]=y1[x];

Drugie rozwiązanie można znaleźć z poniższej formuły.

aa=yp1[x]∫((Exp[-∫p[x]\[DifferentialD]x])/((yp1[x])^2))\[DifferentialD]x//Simplify;
yr2[x_]=Normal@Collect[aa,Log[x]];
y2[x_]=aa;

Sprawdzamy drugie rozwiązanie.

y[x_]=aa;
Operator3[y,x]

Wrońskian

Det[{{y1[x],y2[x]},{y1'[x],y2'[x]}}]//Simplify

Rozwiązania szczególne tworzą układ liniowo niezależny, gdyż Wrońskian jest różny od zera.

Podsumowanie

Print["(Subscript[y, 1])(x) = "yr1[x]]

Print["(Subscript[y, 2])(x) = "yr2[x]]

Zatem otrzymaliśmy dwa liniowo niezależne rozwiązania. 
Zobaczmy jak wyglądają uzyskane rozwiązania.

Manipulate[Plot[Evaluate[{yr1[x],yr2[x]}],{x,a,b},PlotStyle->{Red,Green}],{a,(10^(-6)),10},{b,(10^(-5)),10}]

Literatura

Weber H.J., Arfken G.B., Essential mathematical methods for physicists, Academic Press, Amsterdam 2003

Riley K.S., Hobson M.P., Mathematical methods for physics and engineering, Cambridge University Press, Cambridge 1998

Andrzej Lenda, Bartłomiej Spisak, Wybrane rozdziały matematycznych metod fizyki. Rozwiązane problemy, Uczelniane Wydawnictwa Naukowo-Dydaktyczne AGH, Kraków 2006

Andrzej Lenda, Wybrane rozdziały matematycznych metod fizyki
http://www.ftj.agh.edu.pl/~lenda/mmf/skrypt.pdf

Metoda Frebeniusa na mathworld.wolfram.com
http://mathworld.wolfram.com/FrobeniusMethod.html