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

Z Wikibooks, biblioteki wolnych podręczników.

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

Rozwiązać numerycznie zagadnienie Cauchy'ego:

na przedziale . Zaimplementować schemat Eulera i schemat Mid-Point.

[edytuj] Konstrukcja schematu Eulera

Załóżmy, że . Wówczas ze wzoru Taylora:

x(t + h) = x(t) + hx'(t) + o(h²) = x(t) + hf(x) + O(h²)

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₀ i dostajemy wszystkie następne wyrazy)
• rzędu h¹ (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 1)

[edytuj] Konstrukcja schematu Mid-point

Załóżmy, że . Wówczas podobnie:

x(t + h) = x(t) + hx'(t) + h²x''(t) / 2 + O(h³)

x(t − h) = x(t) − hx'(t) + h²x''(t) / 2 + O(h³)

Odejmujemy stronami i dostajemy:

x(t + h) = x(t − h) + 2hx'(t) + O(h³)

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₀ i x₁ żeby wystartować)
• rzędu h² (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 2)

[edytuj] 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 , gdzie N = (b − a) / h. Zadajemy x₀ = 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₁ = x₀ + hf(x₀). Kolejne wyrazy obliczamy według reguły podanej w schemacie.
3. Rozwiązaniem jest funkcja, przybliżona w punktach: x(a + kh) = x_k

[edytuj] 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;

[edytuj] 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;

[edytuj] 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:

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.

[edytuj] 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":

na odcinku [0,10]. 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 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.