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{array}{c}{x}^{\prime }=f(x)\\ x(a)=x_0\end{array}$

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)+h{x}^{\prime }(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)+h{x}^{\prime }(t)+h^2{x}^{″}(t)/2+O(h^3)$

$x(t-h)=x(t)-h{x}^{\prime }(t)+h^2{x}^{″}(t)/2+O(h^3)$

Odejmujemy stronami i dostajemy:

$x(t+h)=x(t-h)+2h{x}^{\prime }(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

Plik:Schemat Eulera a Midpoint.png

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

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

$\{\begin{array}{c}{x}^{\prime }=f(x)=-x,t\in [0,5]\\ x(0)=1\end{array}$

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{array}{c}{x}^{\prime }=998x+1998y\\ y^{\prime }=-999x-1999y\\ x(0)=y(0)=1\end{array}$

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

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

$A=\left(\begin{array}{cc}998& 1998\\ -999& -1999\end{array}\right)$

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.

Szablon:Nawigacja