Przejdź do zawartości

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

Z Wikibooks, biblioteki wolnych podręczników.

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:

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

  1. , 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 .