Metoda Różnic Skończonych

Równania różniczkowe i całkowe. Równania różnicowe. Transformata Laplace'a i Fouriera oraz ich zastosowanie w równaniach różniczkowych.
janusz47
Użytkownik
Użytkownik
Posty: 7910
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1670 razy

Metoda Różnic Skończonych

Post autor: janusz47 »

Mamy do rozwiązania zagadnienie początkowe równania różniczkowego zwyczajnego drugiego rzędu:

\(\displaystyle{ y''(x) +p(x)y'+q(x)y = r(x), \ \ a \leq x \leq b \ \ (1) }\)

\(\displaystyle{ y(a) = \alpha, \ \ y(b) = \beta.}\)

Dzielimy przedział \(\displaystyle{ [a,\ \ b] }\) na \(\displaystyle{ N }\) równych podprzedziałów, każdy o długości \(\displaystyle{ h, }\)

\(\displaystyle{ h = \frac{b-a}{N}.}\)

Aproksymujemy równanie \(\displaystyle{ (1) }\) różnicami centralnymi:

\(\displaystyle{ y'(x_{i}) \approx \frac{1}{2h}[y(x_{i+1} -y(x_{i-1})] \ \ (2) }\)

\(\displaystyle{ y''(x_{i}) \approx \frac{1}{h^2} [y(x_{i+1} - 2y(x_{i}) + y(x_{i-1})] \ \ (3)}\)

Podstawiając \(\displaystyle{ (2), (3) }\) do \(\displaystyle{ (1) }\) otrzymujemy równanie różnicowe

\(\displaystyle{ \left[ 1 -\frac{h}{2 }p_{i}\right] y_{i-1} + (-2 +h^2 q_{i}) y_{i} + \left[ 1 + \frac{h}{2}p_{i}\right] y_{i+1}= h^2 r_{i} \ \ (4) }\)

\(\displaystyle{ i=1,2,...,N-1 }\) i \(\displaystyle{ y_{0} = \alpha, \ \ y_{N}=\beta, y_{i}\approx y(x_{i}), \ \ p_{i}=p(x_{i}), \ \ q_{i}=q(x_{i}), \ \ r_{i}= r(x_{i}). }\)

Równanie \(\displaystyle{ (4) }\) zapisujemy w postaci macierzowej

\(\displaystyle{ A\vec{y} = \vec{b} }\)

gdzie:

\(\displaystyle{ A = \left [\begin{matrix} b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} & 0 \\ & & \ddots & &\\ & 0 & a_{N-2} & b_{N-2}& c_{N-2}\\ & & & a_{N-1} & b_{N-1} \end{matrix}\right], }\)

\(\displaystyle{ \vec{y} = \left[ \begin{matrix} y_{1}\\y_{2}\\ \vdots \\ y_{N-2}\\ y_{N-1}\end{matrix} \right],}\)

\(\displaystyle{ \vec{b} =\left[\begin{matrix} d_{1}-a_{1}\alpha \\ d_{2} \\ \vdots \\ d_{N-2} \\ d_{N-1}-c_{N-1}\beta \end{matrix} \right] }\)

Wielkości \(\displaystyle{ a_{i}, \ \ b_{i}, \ \ c_{i}, \ \ d_{i} }\) definiujemy równaniami:

\(\displaystyle{ a_{i} = 1 -\frac{h}{2}p_{i}, \ \ b_{i}= -2 + h^2q_{i}, \ \ c_{i}=1 + \frac{h}{2}p_{i}, \ \ d_{i} = h^2r_{i}.}\)


Przykład

Rozwiązanie równania

\(\displaystyle{ y'' +e^{x}y' -xy =e^{-x}(-x^2-2x-3) -x+2 }\)

z warunkami brzegowymi: \(\displaystyle{ y(0)=-1, \ \ y(1) = 0.}\)

Program w MATLAB

Kod: Zaznacz cały

function mrs(p,q,r,aa,bb,y0,yn,n) 
% Rozwiązanie równania różniczkowego rzędu drugiego. 
fprintf('\n') 
h=(bb-aa)/n; 
blank=' '; 
disp(['Metoda Różnic Skończonych h=',num2str(h)]) 
fprintf('\n') 
for i=1:n-1 
 x=aa+i*h; 
 if (i~=1) 
 a(i-1)=1-h/2*feval(p,x); %Obliczenie poddiagonali macierzy A. 
 end 
 b(i)=-2+h^2*feval(q,x); %Obliczenie głównej diagonali macierzy A. 
 if (i~=n-1) 
 c(i)=1+h/2*feval(p,x); %Obliczenie naddiagonali macierzy A. 
 end 
end 
disp(' Poddiagonala macierzy   A =') 
disp(a) 
disp(' Główna diagonala macierzy A =') 
disp(b) 
disp(' Naddiagonala macierzy  A =') 
disp(c) 
% Obliczenie współrzędnych wektora b 
d(1)=h^2*feval(r,aa+h)-y0*(1-h/2*feval(p,aa+h)); 
d(n-1)=h^2*feval(r,bb-h)-yn*(1+h/2*feval(p,bb-h)); 
for i=2:n-2 
 x=aa+i*h; 
 d(i)=h^2*feval(r,x); 
end 
fprintf('\n') 
disp(' Współczynnik b =') 
disp(d) 
disp(' Rozwiązanie  =') 
fprintf('\n') 
disp('   xi           yi     ') 
fprintf('\n') 
% Rozwiązanie układu trójdiagonalnego 
for i=2:n-1 
 ymult=a(i-1)/b(i-1); 
 b(i)=b(i)-ymult*c(i-1); 
 d(i)=d(i)-ymult*d(i-1); 
end 
y(n)=yn; 
y(n-1)=d(n-1)/b(n-1); 
for i=n-2:-1:1 
 y(i)=(d(i)-c(i)*y(i+1))/b(i); 
end 
fprintf('%6.2f %12.6f\n',aa,y0) 
for i=1:n 
 x=aa+i*h; 
  fprintf('%6.2f %12.6f\n',x,y(i)) 
  end 
Przykład

Kod: Zaznacz cały


function p =p1(x)
p=exp(x);

function q=q1(x)
q=-x;

function r =r1(x)
r=exp(-x)*(-x.^2+2*x-3)-x+2;
 mrs('p1','q1','r1',0,1,-1,0,10)

Metoda Różnic Skończonych h= 0.1

 Poddiagonala macierzy   A =
    0.9389    0.9325    0.9254    0.9176    0.9089    0.8993    0.8887    0.8770

 Główna diagonala macierzy A =
   -2.0010   -2.0020   -2.0030   -2.0040   -2.0050   -2.0060   -2.0070   -2.0080   -2.0090

 Naddiagonala macierzy  A =
    1.0553    1.0611    1.0675    1.0746    1.0824    1.0911    1.1007    1.1113


 Współczynnik b =
    0.9383   -0.0036   -0.0014    0.0002    0.0014    0.0021    0.0026    0.0028    0.0028

 Rozwiązanie  =

   xi           yi     

  0.00    -1.000000
  0.10    -0.814197
  0.20    -0.654714
  0.30    -0.518229
  0.40    -0.401815
  0.50    -0.302890
  0.60    -0.219182
  0.70    -0.148691
  0.80    -0.089661
  0.90    -0.040549
  1.00     0.000000

Ostatnio zmieniony 11 sie 2022, o 19:29 przez Jan Kraszewski, łącznie zmieniany 1 raz.
Powód: Poprawa wiadomości.
ODPOWIEDZ