GNU Octave/Numeryczne rozwiązywanie równań różniczkowych
Numeryczne rozwiązywanie równań różniczkowych zwyczajnych pierwszego rzędu
edytujZa 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:
- , gdzie
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ąż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
edytujRozwią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ć:
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
edytujRozpatrzmy 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ę .
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
edytujZnaleźć 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 :
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.
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 .