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

Z Wikibooks, biblioteki wolnych podręczników.
Przejdź do nawigacji Przejdź do wyszukiwania

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

Rozwiązać numerycznie zagadnienie Cauchy'ego:

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

Konstrukcja schematu Eulera[edytuj]

Załóż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[edytuj]

Załóż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[edytuj]

Rozwiązanie zagadnienia Cauchy'ego na przedziale za pomocą pewnego schematu:

  1. Wyznaczamy krok , który powinien być "wystarczająco" mały.
  2. 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.
  3. 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[edytuj]

Porównanie schematów Eulera i Midpoint dla funkcji .

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

Rozwią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.