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

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.