GNU Octave/Rozwiązywanie równań różniczkowych

Rozwiązywanie równań różniczkowych implementując schematy różnicowe

Rozwiązać numerycznie zagadnienie Cauchy'ego:

${\begin{cases}x'=f(x)\\x(a)=x_{0}\end{cases}}$

na przedziale $t\in [a,b]$. Zaimplementować schemat Eulera i schemat Mid-Point.

Konstrukcja schematu Eulera

Załóżmy, że $x\in C^{1}[a,b]$. Wówczas ze wzoru Taylora:

$x(t+h)=x(t)+hx'(t)+o(h^{2})=x(t)+hf(x)+O(h^{2})$

Oznaczamy $x(t+kh)=x_{k}$ i pomijamy wyrazy drugiego rzędu:

$x_{k}=x_{k-1}+hf(x_{k})$

Schemat ten jest:

• otwarty (kolejny wyraz wylicza się jako nieuwikłana funkcja poprzednich)
• jednokrokowy (do obliczenia $x_{k+1}$ używa tylko $x_{k}$)
• samostartujący (podajemy $x_{0}$ i dostajemy wszystkie następne wyrazy)
• rzędu $h^{1}$ (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 1)

Konstrukcja schematu Mid-point

Załóżmy, że $x\in C^{2}[a,b]$. Wówczas podobnie:

$x(t+h)=x(t)+hx'(t)+h^{2}x''(t)/2+O(h^{3})$

$x(t-h)=x(t)-hx'(t)+h^{2}x''(t)/2+O(h^{3})$

Odejmujemy stronami i dostajemy:

$x(t+h)=x(t-h)+2hx'(t)+O(h^{3})$

Pomijamy wyrazy rzędu 3 i dyskretyzujemy:

$x_{k+1}=x_{k-1}+2hf(x_{k})$

Schemat ten jest:

• otwarty
• dwukrokowy (do obliczenia $x_{k+1}$ używa $x_{k}$ i $x_{k-1}$)
• nie samostartujący (potrzebuje $x_{0}$ i $x_{1}$ żeby wystartować)
• rzędu $h^{2}$ (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 $[a,b]$ za pomocą pewnego schematu:

1. Wyznaczamy krok $h$, który powinien być "wystarczająco" mały.
2. Obliczamy ciąg $(x_{k})_{k=0}^{k=N}$, gdzie $N=(b-a)/h$. Zadajemy $x_{0}=x(a)$. Jeśli schemat nie jest samostartujący możemy zadać brakujące wyrazy początkowe wyznaczone na przykład za pomocą schematu Eulera, tj. $x_{1}=x_{0}+hf(x_{0})$. Kolejne wyrazy obliczamy według reguły podanej w schemacie.
3. Rozwiązaniem jest funkcja, przybliżona w punktach: $x(a+kh)=x_{k}$

Implementacja schematu Eulera

 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

 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

Porównanie schematów Eulera i Midpoint dla funkcji $x(t)=exp(-t)$.

Przetestujmy schematy Euler i Mid-point dla równania:

${\begin{cases}x'=f(x)=-x,t\in [0,5]\\x(0)=1\end{cases}}$

Rozwiązaniem tego równania jest $x(t)=\exp(-t)$.

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

Rozwiązać za pomocą funkcji lsode sztywny układ równań używając schematu Adamsa i "stiff":

${\begin{cases}x'=998x+1998y\\y'=-999x-1999y\\x(0)=y(0)=1\end{cases}}$

na odcinku $[0,10]$. Sprawdzić czy ten układ sztywny?

Sztywność układu sprawdzamy obliczając wartości własne macierzy:

$A={\begin{pmatrix}998&1998\\-999&-1999\end{pmatrix}}$

Dostajemy:

 A=[998 1998; -999 -1999];
 m=eig(A)
 m =
     -1
  -1000

a zatem współczynnik sztywności wynosi wynosi $1000$, 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.