Praktická cvičení z Numerické matematiky

Okrajová úloha pro ODR, speciálně Dirichletova okrajová podmínka

  • Dirichletova okrajová úloha, p(x) = 1, q(x) = 0. Cely skript zde: bvp4odr.m
    • Rovnice           - y ' ' = f       v intervalu     I = (a, b)
    Dirichletovy okrajové podmínky: y(a) = Ya, y(b) = Yb 

    n  = 20;
a  = 0; b  = 1;     h = (b - a) / n;
Ya = 0; Yb = 0;  % okrajova podminka

x  = (a:h:b)';
x0 = x(2:end-1);
x2 = (x(1:end-1) + x(2:end)) / 2.;

A = diag(ones(n - 1, 1) * 2);
A = A - diag(ones(n - 2, 1), -1);
A = A - diag(ones(n - 2, 1),  1);

fx = exp(-(x0 - 0.4).^2);
B = h * h * fx; 
B(1)   = B(1)   + Ya;
B(end) = B(end) + Yb;

Y = A \ B;
Y = [Ya; Y; Yb];
  • Dirichletova/Neumannova/Sturmova okrajová úloha, p(x) = 1, q(x) = 0. Cely skript zde: bvp4odr.m
    • Rovnice           - y ' ' = f       v intervalu     I = (0,1)
    Dirichletova okrajová podmínka v x = a: y(a) = Ya
    Sturmova okrajová podmínka v x = b: beta1 y'(b) + beta2 y(b) = beta2 Yb, :

    Aproximace (beta1 + beta2 * h) y(b) - beta1 * y(b - h) = h beta2 Yb: 


    n  = 20;
a  = 0; b  = 1;     h = (b - a) / n;
Ya = 0; Yb = 0;  % okrajova podminka

x  = (a:h:b)';
x0 = x(2:end-1);
x2 = (x(1:end-1) + x(2:end)) / 2.;

A = diag(ones(n - 1, 1) * 2);
A = A - diag(ones(n - 2, 1), -1);
A = A - diag(ones(n - 2, 1),  1);

% Sturmova okrajova podminka:
A(n, n) = beta1 + beta2 * h;
A(n, n - 1) = beta1;


fx = exp(-(x0 - 0.4).^2);
B = h * h * fx; 
B(1)   = B(1)   + Ya;
B(end + 1) = h * beta2 * Yb;

Y = A \ B;
Y = [Ya; Y];
  • Úloha pro obecná p(x), q(x), f(x)
    • Zde uvažujeme úlohu - ( p(x) y ') ' + q(x) y = f,     y(0) = 0,   y(1) = 0 

    n  = 10;
a  = 0; b  = 1;     h = (b - a) / 10;
Ya = 0; Yb = -0.1;  % okrajova podminka

x  = (a:h:b)';
x0 = x(2:end-1);
x2 = (x(1:end-1) + x(2:end)) / 2.;


px2 = 1 + x2.^2;             % p(x)
qx  = 0 * x0; %0.1 + 0 * x0.^2;       % q(x)
fx  = sin(x0 * pi);          % f(x)


A = diag(px2(1:end-1) + px2(2:end) + h * h * qx);
A = A - diag(px2(2:end-1), -1);
A = A - diag(px2(2:end-1),  1);

B = h * h * fx; 
B(1)   = B(1)   + px2(1)   * Ya;
B(end) = B(end) + px2(end) * Yb;

Y = A \ B;
Y = [Ya; Y; Yb];