Tekst źródłowy strony GNU Octave/Rozwiązywanie równań różniczkowych

===Rozwiązywanie równań różniczkowych implementując schematy różnicowe===
Rozwiązać numerycznie [[w:zagadnienie Cauchy'ego|zagadnienie Cauchy'ego]]:
:<math>\begin{cases}
x'=f(x)\\
x(a)=x_0
\end{cases}</math>
na przedziale <math>t\in[a,b]</math>. Zaimplementować schemat Eulera i schemat Mid-Point.

====Konstrukcja schematu Eulera====
Załóżmy, że <math>x \in C^1[a,b]</math>. Wówczas ze [[w:wzór Taylora|wzoru Taylora]]:
:<math>x(t+h) = x(t) + h x'(t) + o(h^2) = x(t) + h f(x) + O(h^2) </math>
Oznaczamy <math>x(t+kh)=x_k</math> i pomijamy wyrazy drugiego rzędu:
:<math>x_k=x_{k-1}+h f(x_k)</math>
Schemat ten jest:
*''otwarty'' (kolejny wyraz wylicza się jako nieuwikłana funkcja poprzednich)
*''jednokrokowy'' (do obliczenia <math>x_{k+1}</math> używa tylko <math>x_{k}</math>)
*''samostartujący'' (podajemy <math>x_0</math> i dostajemy wszystkie następne wyrazy)
*''rzędu <math>h^1</math>'' (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 1)

====Konstrukcja schematu Mid-point====
Załóżmy, że <math>x \in C^2[a,b]</math>. Wówczas podobnie:
:<math>x(t+h) = x(t) + h x'(t) + h^2 x''(t) / 2 + O(h^3)</math>
:<math>x(t-h) = x(t) - h x'(t) + h^2 x''(t) / 2 + O(h^3)</math>
Odejmujemy stronami i dostajemy:
:<math>x(t+h)=x(t-h)+2 h x'(t) + O(h^3)</math>
Pomijamy wyrazy rzędu 3 i dyskretyzujemy:
:<math>x_{k+1}=x_{k-1}+2 h f(x_k)</math>
Schemat ten jest:
*''otwarty''
*''dwukrokowy'' (do obliczenia <math>x_{k+1}</math> używa <math>x_k</math> i <math>x_{k-1}</math>)
*''nie samostartujący'' (potrzebuje <math>x_0</math> i <math>x_1</math> żeby wystartować)
*''rzędu <math>h^2</math>'' (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 2)

====Rozwiązywanie równania za pomocą schematu====
Rozwiązanie zagadnienia Cauchy'ego na przedziale <math>[a, b]</math> za pomocą pewnego schematu:
# Wyznaczamy krok <math>h</math>, który powinien być "wystarczająco" mały.
# Obliczamy ciąg <math>(x_k)_{k=0}^{k=N}</math>, gdzie <math>N=(b-a)/h</math>. Zadajemy <math>x_0=x(a)</math>. Jeśli schemat nie jest samostartujący możemy zadać brakujące wyrazy początkowe wyznaczone na przykład za pomocą schematu Eulera, tj. <math>x_1=x_0+h f(x_0)</math>. Kolejne wyrazy obliczamy według reguły podanej w schemacie.
# Rozwiązaniem jest funkcja, przybliżona w punktach: <math>x(a+k h)=x_k</math>

====Implementacja schematu Eulera====
<syntaxhighlight lang="matlab">
 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;
</syntaxhighlight>

====Implementacja schematu Mid-point====
<syntaxhighlight lang="matlab">
 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;
</syntaxhighlight>

====Testowanie schematów====
[[Grafika:Schemat Eulera a Midpoint.png|thumb|right|400px|Porównanie schematów Eulera i Midpoint dla funkcji <math>x(t)=exp(-t)</math>.]]
Przetestujmy schematy Euler i Mid-point dla równania:
:<math>\begin{cases}
x'=f(x)=-x, t \in [0, 5]\\
x(0)=1
\end{cases}</math>
Rozwiązaniem tego równania jest <math>x(t)=\exp(-t)</math>.

Zdefiniujmy funkcję:
<syntaxhighlight lang="matlab">
 function y=fun(x)
    y=-x; 
 endfunction;
</syntaxhighlight>
Zadajmy schematy:
<syntaxhighlight lang="matlab">
 a=0
 b=5
 x0=1
 N=100
</syntaxhighlight>
Obliczmy schematy:
<syntaxhighlight lang="matlab">
 z=midpoint ("fun",a,b,x0,N);
 y=euler    ("fun",a,b,x0,N);
</syntaxhighlight>
Obliczmy rozwiązanie teoretyczne:
<syntaxhighlight lang="matlab">
 t=linspace (a,b,N);
 w=1./exp(t);
</syntaxhighlight>
Rysujmy rozwiązanie teoretyczne i dwa rozwiązania obliczone za pomocą schematów:
<syntaxhighlight lang="matlab"> plot(t,z,"-r;midpoint;",t,y,"-b;euler;",t,w,"-g;rozwiazanie;"); </syntaxhighlight>
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===
Rozwiązać za pomocą funkcji '''lsode''' [[w:Sztywny układ równań|sztywny układ równań]] używając [[w:Schemat Adamsa|schematu Adamsa]] i "stiff":
:<math>\begin{cases}
x' =  998 x + 1998 y \\
y' = -999 x - 1999 y \\
x(0) = y(0) = 1 
\end{cases}</math>
na odcinku <math>[0,10]</math>. Sprawdzić czy ten układ sztywny?

Sztywność układu sprawdzamy obliczając wartości własne macierzy:
:<math>A=\begin{pmatrix}
998 & 1998 \\
-999 &-1999
\end{pmatrix}</math>
Dostajemy:
<syntaxhighlight lang="matlab">
 A=[998 1998; -999 -1999];
 m=eig(A)
 m =
     -1
  -1000
</syntaxhighlight>
a zatem współczynnik sztywności wynosi wynosi <math>1000</math>, czyli możemy określić układ jako sztywny.

Rozwiążmy go więc napierw metodą "stiff". Zdefiniujmy odpowiednią funkcję:
<syntaxhighlight lang="matlab">
 function y=fun_sztywna(x)
    y(1)= 998*x(1)+1998*x(2);
    y(2)=-999*x(1)-1999*x(2);
 endfunction;
</syntaxhighlight>
Rozwiązujemy i mierzymy czas:
<syntaxhighlight lang="matlab">
 lsode_options("integration method", "stiff");
 tic; z=lsode("fun_sztywna", x0, t); toc
 ans = 0.32850
</syntaxhighlight>
Spróbujmy także schematem Adamsa:
<syntaxhighlight lang="matlab">
 lsode_options("integration method", "adams");
 tic; z=lsode("fun_sztywna", x0, t); toc
 ans = 26.846
</syntaxhighlight>

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.

<noinclude>
{{Nawigacja|GNU Octave|
[[../Numeryczne rozwiązywanie równań różniczkowych/]]|
[[../Chaos w atraktorze Lorentza/]]|
}}</noinclude>