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

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

Rozwiązanie równania $y^{\prime }=-ax,y(0)=1,a=0.1$, na przedziale [0,10]

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

1. $\{\begin{array}{c}{y}^{\prime }=-ay,x\in [0,1],a=\pm 0.1,\pm 1,\pm 100\\ y(0)=1\end{array}$
2. $\{\begin{array}{c}{y}^{\prime }=y^2,x\in [0,3]\\ y(0)=1\end{array}$
3. $\{\begin{array}{c}{y}^{\prime }=y\cdot x\sin (ax),x\in [0,20]\\ y(0)=1\end{array}$

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

1. $y=e^{ax}$
2. $y=\frac{1}{1-x}$
3. $y=e^{(\sin (ax)-ax\cos (ax))/a^2}$, gdzie $a=100$

Rozwiązanie równania $y^{\prime }=y^2,y(0)=1$, 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 $y^{\prime }=y\cdot x\sin (x)$ 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 $1/(1-x)$ wybucha w $x=1$.

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

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

$\{\begin{array}{c}{x}^{″}=f(x)=-x\\ x(0)=0\\ x^{\prime }(0)=1\end{array}$

Równanie to można scałkować symbolicznie, a jego rozwiązaniem jest $x(t)=sin(t)$.

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 $\bar{x}=[x_1,x_2]=[x(t),x^{\prime }(t)]$ ma postać:

$\{\begin{array}{c}{x_1}^{\prime }=x_2\\ {x_2}^{\prime }=f(x_1)=-x_1\\ x_1(0)=0\\ x_2(0)=1\end{array}$

Wykres rozwiązania $[x,x^{\prime }]$ równania $x^{″}=-x$ z warunkami początkowymi $x(0)=0$, $x^{\prime }(0)=1$

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 $[1;0]$. Narysujmy wykres $[x_1,x_2]$:

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

Wynikiem jest okrąg, czyli zbiór $\{[\sin (t),\cos (t)]:t\in [0,2\pi ]\}$, co zgadza się z rozwiązaniem teoretycznym.

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

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

$\{\begin{array}{c}{x}^{\prime }=f(x,t),x\in [0,1]\\ x(0)=s\end{array}$

Równanie jest sparametryzowane parametrem $s\in [0,1]$ i dla zadanej funkcji $f$. Oznaczmy $y_s(t)$ rozwiązanie równania dla ustalonego $s$. Dla jakiego $s$ dostaniemy $y_s(1)=\frac{1}{2}$? Innymi słowy, jeśli oznaczymy $F(s)=y_s(1)$ należy rozwiązać równanie:

$F(s)=\frac{1}{2}$

Jak wygląda potok $y_s(t)$ i wykres $F(s)$? Rozważmy funkcję $f(y,t)=sin(y^2+1)*y*(1-y)$.

Wykres potoku dla zagadnienia Cauchy'ego dla równania różniczkowego $x^{\prime }=\sin (x^2+1)x(1-x)$ z warunkiem początkowym $x(0)=s$ na przedziale $t\in [0,1]$ parametryzowanego parametrem $s\in [0,1]$. Linia, która "kończy" się w $(1,0.5)$ zaczyna się w przybliżeniu w $(0,0.3)$.

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

Znaleźć funkcję $u=u(t)$ spełniającą równanie różniczkowe z warunkami brzegowymi:

$\{\begin{array}{c}{u}^{″}+u^2=f(t)\\ u(0)=0,u(1)=0\end{array}$

Dla funkcji $f(t)=sin(t)$.

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

$\{\begin{array}{c}{u}^{″}+u^2=f(t)\\ u(0)=0,u^{\prime }(0)=s_0\end{array}$

spełniało zadany warunek brzegowy w drugim końcu przedziału, tj. $u(1)=0$. Można to zrobić metodą z "odwróconym" warunkiem początkowym, opisaną wcześniej. Przekształćmy zagadnienie na równanie pierwszego rzędu zmiennej $[x_1,x_2]=[u(t),u^{\prime }(t)]$:

$\{\begin{array}{c}{x_1}^{\prime }=x_2\\ {x_2}^{\prime }=-x_1^2+f(t)\\ x_1(0)=0\\ x_2(0)=s_0\end{array}$

Potok fazowy dla równania $u^{″}+u^2=sin(t)$ z warunkami początkowymi $u(0)=0$ i $u^{\prime }(0)=s$ dla parametrów $s\in [-1,1]$.

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 $s_0\in [-1,1]$. (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 $(0,0)$ oraz $(1,0)$ ma w $t=0$ nachylenie bliskie 0.

Rozwiązanie równania $u^{″}+u^2=sin(t)$ z warunkami brzegowymi $u(0)=0$ i $u(1)=0$.

Wyznaczymy teraz tę wartość dokładnie. Zdefiniujmy funkcję $u_s(1)$

 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 $(0,0)$ i $(1,0)$.