GNU Octave: Różnice pomiędzy wersjami
Usunięta treść Dodana treść
schematy różnicowe |
|||
Linia 433:
plot(t, z(:,1);
Rozwiązanie rzeczywiście przechodzi przez punkty <math>(0,0)</math> i <math>(1,0)</math>.
===Rozwiązywanie równań różniczkowych za pomocą schematów różnicowych===
Rozwiązać numerycznie [[zagadnienie Cauchy'ego]]:
:<math>\begin{cases}
x'=f(x)\\
x(a)=x_0
\end{cases}</math>
na przedziale <math>t\in[a,b]</math>, za pomocą schamtu Eulera i schematu Mid-Point.
====Konstrukcja schematu Eulera====
Załóżmy, że <math>x \in C^1[a,b]</math>. Wówczas ze [[wzór Taylora|wzoru Taylora]]:
:<math>x(t+h) = x(t) + h x'(t) + o(h^2) = x(t) + h f(x) + O(h^2) </math>
Oznaczamy <math>x(t+kh)=x_k</math> i pomijamy wyrazy drugiego rzędu:
:<math>x_k=x_{k-1}+h f(x_k)</math>
Schemat ten jest:
*''otwarty'' (kolejny wyraz wylicza się jako nieuwikłana funkcja poprzednich)
*''jednokrokowy'' (do obliczenia <math>x_{k+1}</math> używa tylko <math>x_{k}</math>)
*''samostartujący'' (podajemy <math>x_0</math> i dostajemy wszystkie następne wyrazy)
*''rzędu <math>h^1</math>'' (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 1)
====Konstrukcja schematu Mid-point====
Załóżmy, że <math>x \in C^2[a,b]</math>. Wówczas podobnie:
:<math>x(t+h) = x(t) + h x'(t) + h^2 x''(t) / 2 + O(h^3)</math>
:<math>x(t-h) = x(t) - h x'(t) + h^2 x''(t) / 2 + O(h^3)</math>
Odejmujemy stronami i dostajemy:
:<math>x(t+h)=x(t-h)+2 h x'(t) + O(h^3)</math>
Pomijamy wyrazy rzędu 3 i dyskretyzujemy:
:<math>x_{k+1}=x_{k-1}+2 h f(x_k)</math>
Schemat ten jest:
*''otwarty''
*''dwukrokowy'' (do obliczenia <math>x_{k+1}</math> używa <math>x_k</math> i <math>x_{k-1}</math>)
*''nie samostartujący'' (potrzebuje <math>x_0</math> i <math>x_1</math> żeby wystartować)
*''rzędu <math>h^2</math>'' (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 2)
====Rozwiązywanie równania za pomocą schematu====
Rozwiązanie zagadnienia Cauchy'ego na przedziale <math>[a, b]</math> za pomocą pewnego schematu:
# Wyznaczamy krok <math>h</math>, który powinien być "wystarczająco" mały.
# Obliczamy ciąg <math>(x_k)_{k=0}^{k=N}</math>, gdzie <math>N=(b-a)/h</math>.
Zadajemy <math>x_0=x(a)</math>. Jeśli schemat nie jest samostartujący możemy zadać brakujące
wyrazy początkowe wyznaczone na przykład za pomocą schematu Eulera, tj. <math>x_1=x_0+h f(x_0)</math>.
Kolejne wyrazy obliczamy według reguły podanej w schemacie.
# Rozwiązaniem jest funkcja, przybliżona w punktach:
:<math>x(a+k h)=x_k</math>
====Implementacja schematu Eulera====
function y=euler(fun,a,b,x0,N)
h=(b-a)/(N-1);
y=zeros(N,1);
y(1)=x0;
for (k=2:N)
y(k)=y(k-1)+h*feval(fun,y(k-1));
endfor;
endfunction;
====Implementacja schematu Mid-point====
function y=midpoint(fun,a,b,x0,N)
h=(b-a)/(N-1);
y=zeros(N,1);
y(1)=x0;
y(2)=x0+h*feval(fun,x0);
for (k=3:N)
y(k)=y(k-2)+2*h*feval(fun,y(k-1));
endfor;
endfunction;
====Testowanie schematów====
[[Grafika:Schemat Eulera a Midpoint.png|thumb|right|400px|Porównanie schematów Eulera i Midpoint dla funkcji <math>x(t)=exp(-t)</math>.]]
Przetestujmy schematy Euler i Mid-point dla równania:
:<math>\begin{cases}
x'=f(x)=-x, t \in [0, 5]\\
x(0)=1
\end{cases}</math>
Rozwiązaniem tego równania jest <math>x(t)=\exp(-t)</math>.
Zdefiniujmy funkcję:
function y=fun(x)
y=-x;
endfunction;
Zadajmy schematy:
a=0
b=5
x0=1
N=100
Obliczmy schematy:
z=midpoint ("fun",a,b,x0,N);
y=euler ("fun",a,b,x0,N);
Obliczmy rozwiązanie teoretyczne:
t=linspace (a,b,N);
w=1./exp(t);
Rysujmy rozwiązanie teoretyczne i dwa rozwiązania obliczone za pomocą schematów:
plot(t,z,"-r;midpoint;",t,y,"-b;euler;",t,w,"-g;rozwiazanie;");
Widać, że schemat midpoint lepiej przybliża rozwiązanie na początku (jest wyższego rzędu), ale potem "rozjeżdza" się,
w przeciwieństwie do schematu Eulera, który okazuje się być stabilniejszy.
===Wartości własne===
| |||