\(\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
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