GNU Octave/Rozwiązywanie równań różniczkowych
Rozwiązywanie równań różniczkowych implementując schematy różnicowe
edytujRozwiązać numerycznie zagadnienie Cauchy'ego:
na przedziale . Zaimplementować schemat Eulera i schemat Mid-Point.
Konstrukcja schematu Eulera
edytujZałóżmy, że . Wówczas ze wzoru Taylora:
Oznaczamy i pomijamy wyrazy drugiego rzędu:
Schemat ten jest:
- otwarty (kolejny wyraz wylicza się jako nieuwikłana funkcja poprzednich)
- jednokrokowy (do obliczenia używa tylko )
- samostartujący (podajemy i dostajemy wszystkie następne wyrazy)
- rzędu (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 1)
Konstrukcja schematu Mid-point
edytujZałóżmy, że . Wówczas podobnie:
Odejmujemy stronami i dostajemy:
Pomijamy wyrazy rzędu 3 i dyskretyzujemy:
Schemat ten jest:
- otwarty
- dwukrokowy (do obliczenia używa i )
- nie samostartujący (potrzebuje i żeby wystartować)
- rzędu (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 2)
Rozwiązywanie równania za pomocą schematu
edytujRozwiązanie zagadnienia Cauchy'ego na przedziale za pomocą pewnego schematu:
- Wyznaczamy krok , który powinien być "wystarczająco" mały.
- Obliczamy ciąg , gdzie . Zadajemy . Jeśli schemat nie jest samostartujący możemy zadać brakujące wyrazy początkowe wyznaczone na przykład za pomocą schematu Eulera, tj. . Kolejne wyrazy obliczamy według reguły podanej w schemacie.
- Rozwiązaniem jest funkcja, przybliżona w punktach:
Implementacja schematu Eulera
edytuj function y=euler(fun,a,b,x0,N)
h=(b-a)/(N-1);
y=zeros(N,1);
y(1)=x0;
for (k=2:N)
y(k)=y(k-1)+h*feval(fun,y(k-1));
endfor;
endfunction;
Implementacja schematu Mid-point
edytuj function y=midpoint(fun,a,b,x0,N)
h=(b-a)/(N-1);
y=zeros(N,1);
y(1)=x0;
y(2)=x0+h*feval(fun,x0);
for (k=3:N)
y(k)=y(k-2)+2*h*feval(fun,y(k-1));
endfor;
endfunction;
Testowanie schematów
edytujPrzetestujmy schematy Euler i Mid-point dla równania:
Rozwiązaniem tego równania jest .
Zdefiniujmy funkcję:
function y=fun(x)
y=-x;
endfunction;
Zadajmy schematy:
a=0
b=5
x0=1
N=100
Obliczmy schematy:
z=midpoint ("fun",a,b,x0,N);
y=euler ("fun",a,b,x0,N);
Obliczmy rozwiązanie teoretyczne:
t=linspace (a,b,N);
w=1./exp(t);
Rysujmy rozwiązanie teoretyczne i dwa rozwiązania obliczone za pomocą schematów:
plot(t,z,"-r;midpoint;",t,y,"-b;euler;",t,w,"-g;rozwiazanie;");
Widać, że schemat midpoint lepiej przybliża rozwiązanie na początku (jest wyższego rzędu), ale potem "rozjeżdza" się, w przeciwieństwie do schematu Eulera, który okazuje się być stabilniejszy.
Rozwiązywanie równań różniczkowych z użyciem gotowych schematów różnicowych
edytujRozwiązać za pomocą funkcji lsode sztywny układ równań używając schematu Adamsa i "stiff":
na odcinku . Sprawdzić czy ten układ sztywny?
Sztywność układu sprawdzamy obliczając wartości własne macierzy:
Dostajemy:
A=[998 1998; -999 -1999];
m=eig(A)
m =
-1
-1000
a zatem współczynnik sztywności wynosi wynosi , czyli możemy określić układ jako sztywny.
Rozwiążmy go więc napierw metodą "stiff". Zdefiniujmy odpowiednią funkcję:
function y=fun_sztywna(x)
y(1)= 998*x(1)+1998*x(2);
y(2)=-999*x(1)-1999*x(2);
endfunction;
Rozwiązujemy i mierzymy czas:
lsode_options("integration method", "stiff");
tic; z=lsode("fun_sztywna", x0, t); toc
ans = 0.32850
Spróbujmy także schematem Adamsa:
lsode_options("integration method", "adams");
tic; z=lsode("fun_sztywna", x0, t); toc
ans = 26.846
Jak widać mimo, że schemat "stiff" jest schematem zamkniętym (w każdym kroku rozwiązuje równanie), to działa on znacznie szybciej niż otwarty schemat Adamsa, gdyż rekompensuje sobie długością kroku.