GNU Octave/Numeryczne rozwiązywanie równań różniczkowych

Numeryczne rozwiązywanie równań różniczkowych zwyczajnych pierwszego rzędu

edytuj
 
Rozwiązanie równania  , na przedziale [0,10]

Za pomocą funkcji lsode rozwiązać równania różniczkowe, tj. znaleźć funkcję y=y(x), taką, że:

  1.  
  2.  
  3.  

Rozwiązania powyższych równań można wyznaczyć symbolicznie, są to:

  1.  
  2.  
  3.  , gdzie  
 
Rozwiązanie równania  , na przedziale [0, 0.99]

Zdefiniujmy funkcje, które występują w równaniach:

 function y=linfun(x)
    global a;
    y=a.*x;
 endfunction;

 function y=kwadfun(x)
    y=x.*x;
 endfunction;

 function y=mixfun(x,t)
    y=x.*t.*sin(t*100);
 endfunction;
 
Rozwiązanie równania   na przedziale [0,20]. Dosalismy fale o bardzo dużym nachyleniu i coraz większych amplitudach.

Rozwiążmy pierwsze równanie: Zadajemy wartość dla współczynnika a:

 global a;
 a=0.1;
 %Zadajemy przedział, na którym będziemy rozwiązywać równanie:
  t=[0:0.01:1];
 %Zadajemy warunek brzegowy w pierwszym punkcie czasowym <tt>t(1)</tt>, czyli w <tt>t=0</tt>:
  x0=1; 
 %Rozwiązujemy równanie:
  Z=lsode("linfun", x0, t);
 %Rysujemy rozwiązanie:
  plot(t,Z,"-r")

Podobnie rozwiązujemy drugie równanie:

 t=[0:0.01:0.99];
 x0=1;
 Z=lsode("kwadfun", x0, t);

Funkcja   wybucha w  .

Rozwiązujemy trzecie równanie:

 t=[0:0.01:20];
 x0=1;
 Z=lsode("mixfun", x0, t);
 axis([0, 20, 0, 2]);
 plot(t, Z, "-r");

Numeryczne rozwiązywanie równań różniczkowych zwyczajnych wyższych rzędów

edytuj

Rozwiązać równanie różniczkowe 2. rzędu:

 

Równanie to można scałkować symbolicznie, a jego rozwiązaniem jest  .

Aby rozwiązać równanie numerycznie, przekształcimy je do równania pierwszego rzędu, w którym zmienna jest wektorem i rozwiążemy za pomocą funkcji lsode. Przekształcone równanie dla zmiennej   ma postać:

 
 
Wykres rozwiązania   równania   z warunkami początkowymi  ,  

Zdefiniujmy funkcję:

 function y=ujfun(x)
    y(1)= x(2);
    y(2)=-x(1);
 endfunction;
%Rozwiązujemy równanie <math>\overline{x}'=\texttt{ujfun}(\overline{x})</math>: 
 x0=[1;0];
 t=[0:0.005:100];
 Z=lsode("ujfun", x0, t);

Zauważmy, że warunek początkowy jest pionowym 2-elementowym pionowym wektorem  . Narysujmy wykres  :

 plot(Z(:,2),Z(:,1),"-r");

Wynikiem jest okrąg, czyli zbiór  , co zgadza się z rozwiązaniem teoretycznym.

Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z "odwróconym" warunkiem początkowym

edytuj

Rozpatrzmy zagadnienie Cauchy'ego dla równania różniczkowego:

 

Równanie jest sparametryzowane parametrem   i dla zadanej funkcji  . Oznaczmy   rozwiązanie równania dla ustalonego  . Dla jakiego   dostaniemy  ? Innymi słowy, jeśli oznaczymy   należy rozwiązać równanie:

 

Jak wygląda potok   i wykres  ? Rozważmy funkcję  .

 
Wykres potoku dla zagadnienia Cauchy'ego dla równania różniczkowego   z warunkiem początkowym   na przedziale   parametryzowanego parametrem  . Linia, która "kończy" się w   zaczyna się w przybliżeniu w  .

Najpierw rozwiążemy szukane równanie. Zdefiniujmy funkcję występującą w równaniu:

 function y=dfun(x)
    y=sin(x.^2.+1).*x.*(1.0.-x);
 endfunction;

%Rozwiązujemy równanie dla zadanego parametru <math>s</math> na przedziale <math>[0,1]</math>:
 function y=dfunSolveODE(s)
    t=[0:0.01:1];
    y=lsode("dfun", s, t);
 endfunction;

%Rysujemy wykres funkcji <math>F(s)=y_s(1)=dfunsolve(s)</math>:
 s=[0:0.01:1];
 y=dfunSolveODE(s);
 plot(s,y);

%Rozwiązujemy równanie <math>F(s)-\frac{1}{2}=0</math>. Zdefiniujmy funkcję:
 function y=dfunMinusPol(s)
    z=dfunSolveODE(s);
    y=z(length(z))-0.5;
 endfunction;

% i znajdźmy jej [[w:miejsce zerowe|miejsce zerowe]]:
 fsolve("dfunMinusPol",0.3)
 ans = 0.28628

Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z warunkami brzegowymi

edytuj

Znaleźć funkcję   spełniającą równanie różniczkowe z warunkami brzegowymi:

 

Dla funkcji  .

Aby rozwiązać to zagadnienie trzeba znaleźć taki parametr  , żeby rozwiązanie zagadnienia

 

spełniało zadany warunek brzegowy w drugim końcu przedziału, tj.  . Można to zrobić metodą z "odwróconym" warunkiem początkowym, opisaną wcześniej. Przekształćmy zagadnienie na równanie pierwszego rzędu zmiennej  :

 
 
Potok fazowy dla równania   z warunkami początkowymi   i   dla parametrów  .

Zdefiniujmy funkcje występujące w równaniu:

 function y=d1fun(x,t)
    y(1)=x(2);
    y(2)=-x(1)^2+d1funF(t);
 endfunction;
 function y=d1funF(t)
    y=sin(t);
 endfunction;

%Zdefiniujmy funkcję, która rozwiązuje równanie dla zadanego parametru:
 function y=d1funSolveODE(s)
    t=[0:0.01:1];
    x0=[0;s];
    y=lsode("d1fun", x0, t);
 endfunction;

Narysujmy potok fazowy dla tego równania dla wartości  . (Dokładniej: Funkcja d1funSolveODE działa poprawnie dla pojedynczej wartości parametru s. Aby wywołać rysowanie tak, jak poniżej, trzeba lekko zmienić jej treść).

 s=[-1:0.03:1];
 y=d1funSolveODE(s);
 plot(t,y);

Widzimy, że szukana krzywa przechodząca przez punkty   oraz   ma w   nachylenie bliskie 0.

 
Rozwiązanie równania   z warunkami brzegowymi   i  .

Wyznaczymy teraz tę wartość dokładnie. Zdefiniujmy funkcję  

 function y=d1funFunIn1(s)
    z=d1funSolveODE(s);
    w=z(:,1);
    y=w(length(w));
 endfunction;
% i rozwiążmy równanie:
 s0=fsolve("d1funFunIn1", 0)

%Dostajemy wartość: 
 s0 = -0.15769

%Rysujemy rozwiązanie:
 z=d1funSolveODE(s0);
 plot(t, z(:,1);

Rozwiązanie rzeczywiście przechodzi przez punkty   i  .