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

Z Wikibooks, biblioteki wolnych podręczników.

< GNU Octave

Skocz do: nawigacji, wyszukiwania

Spis treści

1 Numeryczne rozwiązywanie równań różniczkowych zwyczajnych pierwszego rzędu
2 Numeryczne rozwiązywanie równań różniczkowych zwyczajnych wyższych rzędów
3 Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z "odwróconym" warunkiem początkowym
4 Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z warunkami brzegowymi

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

Rozwiązanie równania $y'=-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:

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

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

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

Rozwiązanie równania $y'=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'=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");

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

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

$\begin{cases} x''=f(x)=-x \\ x(0)=0 \\ x'(0)=1 \\ \end{cases}$

Równanie to można scałkować symbolicznie, a jego rozwiązaniem jest $x (t) = s i n (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 $\overline{x}=[x_1, x_2]=[x(t), x'(t)]$ ma postać:

$\begin{cases} x_1' = x_2 \\ x_2' = f(x_1) = -x_1 \\ x_1(0) = 0 \\ x_2(0) = 1 \end{cases}$

Wykres rozwiązania

[x, x']

równania

x'' = - x

z warunkami początkowymi

x (0) = 0

,

x'(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.

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

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

$\begin{cases} x' = f(x,t), x \in [0,1] \\ x(0) = s \end{cases}$

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) = s i n (y 2 + 1) * y * (1 - y)$ .

Wykres potoku dla zagadnienia Cauchy'ego dla równania różniczkowego

x' = 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

[edytuj] 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{cases} u''+u^2=f(t) \\ u(0)=0, u(1)=0 \end{cases}$

Dla funkcji $f (t) = s i n (t)$ .

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

$\begin{cases} u''+u^2=f(t) \\ u(0)=0, u'(0)=s_0 \end{cases}$

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'(t)]$ :

$\begin{cases} x_1' = x_2 \\ x_2' = -x_1^2 + f(t) \\ x_1(0)=0 \\ x_2(0)=s_0 \end{cases}$

Potok fazowy dla równania

u'' + u 2 = s i n (t)

z warunkami początkowymi

u (0) = 0

i

u'(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 = s i n (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)$ .

« Całkowanie numeryczne

Spis treści

Rozwiązywanie równań różniczkowych »

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

Z Wikibooks, biblioteki wolnych podręczników.

Spis treści

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

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

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

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

Widok

Osobiste

Nawigacja

Drukuj lub eksportuj

Szukaj

Narzędzia