Metoda sítí

Lineární ODR 2. řádu (samoadjungovaný tvar)

Předpokládáme úlohu ve tvaru -(p(x) y'(x))' + q(x) y(x) = f(x) na intervalu (a,b), y(a) = y0, y(b) = yn .

Příklad:   -(x² y'(x))' + x³/(2+x) y(x) = x²/(3-x) ,    y(-5) = -2 ,   y(-3) = 0 .

Řešení tohoto příkladu v Matlabu metodou sítí s krokem h = 0.4:

 % sestaveni matice soustavy A
 h = 0.4;                           % volba kroku
 a = -5; b = -3;                    % volba intervalu reseni
 n = round((b-a)/h)-1;              % pocet neznamych
 xh = linspace(a+h, b-h, n);        % vnitrni uzly
 xp = (a + 0.5*h) : h : b;          % deleni intervalu "mezi uzly"
 hq = h^2 * xh .^ 3 ./ (2 + xh);    % hodnoty q v uzlech, vynasobene h^2
 p = xp .^ 2;                       % hodnoty p mezi uzly
 D = hq + p(1:n) + p(2:n+1);        % prvky na hlavni diagonale
 V = - p(2:n);                      % prvky na vedlejsich diagonalach
 A = diag(V,-1) + diag(D) + diag(V,1);  % matice soustavy

 % sestaveni prave strany
 y0 = -2; yn = 0;                % Dirichletovy okrajove podminky
 F = h^2 * xh .^ 2 ./ (3 - xh);  % vypocet funkce f v uzlech, vynasobene h^2
 F(1)=F(1) + p(1)*y0;            % zahrnuti leve okrajove podminky
 F(n)=F(n) + p(n+1)*yn;          % zahrnuti prave okrajove podminky

 % vypocet reseni a nakresleni grafu
 U = A\F';                   % reseni soustavy
 y = [y0; U ; yn];           % pridani okrajovych podminek
 x = [a,xh,b];               % pridani krajnich bodu intervalu
 plot(x,y,x,y,'ro')          % graf numerickeho reseni

Příklad - tlumený oscilátor:   y''(x) + 2 y'(x) + y(x) = e^-x,    y(0) = 2 ,   y(2) = 0 .

Nejprve převedeme úlohu na samoadjungovaný tvar:
-(e^2x y'(x))' - e^2x y(x) = -e^x

Řešení v Matlabu metodou sítí s krokem h = 0.2:

 % sestaveni matice soustavy A
 h = 0.2;                           % volba kroku
 a = 0; b = 2;                      % volba intervalu reseni
 n = round((b-a)/h)-1;              % pocet neznamych
 xh = linspace(a+h, b-h, n);        % vnitrni uzly
 xp = (a + 0.5*h) : h : b;          % deleni intervalu "mezi uzly"
 hq = -h^2 *exp(2*xh);              % hodnoty q v uzlech, vynasobene h^2
 p = exp(2*xp);                     % hodnoty p mezi uzly
 D = hq + p(1:n) + p(2:n+1);        % prvky na hlavni diagonale
 V = - p(2:n);                      % prvky na vedlejsich diagonalach
 A = diag(V,-1) + diag(D) + diag(V,1);  % matice soustavy

 % sestaveni prave strany
 y0 = -2; yn = 0;                % Dirichletovy okrajove podminky
 F = h^2 * xh .^ 2 ./ (3 - xh);  % vypocet funkce f v uzlech, vynasobene h^2
 F(1)=F(1) + p(1)*y0;            % zahrnuti leve okrajove podminky
 F(n)=F(n) + p(n+1)*yn;          % zahrnuti prave okrajove podminky

 % vypocet reseni a nakresleni grafu
 U = A\F';                   % reseni soustavy
 y = [y0; U ; yn];           % pridani okrajovych podminek
 x = [a,xh,b];               % pridani krajnich bodu intervalu
 plot(x,y,'ro')          % graf numerickeho reseni

porovnání výsledků s přesným řešením y = (2 - 2x + 0.5 x²) e^-x

graficky - černá křivka

 xx = linspace(a,b,30);
 yy = (2 - 2*xx + 0.5*xx.^2).*exp(-xx);
 hold on
 plot(xx,yy,'k')

numericky - výpočet maximální chyby

 yy = (2 - 2*x + 0.5*x.^2).*exp(-x);
 max_err = max(abs(y-yy'))