Kvôli tomu, že sme výpočet obmedzili na ohraničený interval, musíme k rovnici dodať okrajové podmienky v krajných bodoch. Krajné body zodpovedajú cenám akcie, ktoré sú veľmi malé, blízke nule (x=-L) a cenám, ktoré sú veľmi veľké a približujú sa k nekonečnu (x=L). Pre takéto limitné hodnoty použijeme aproximácie:
Prvý prístup (explicitná schéma) vyzerá byť jednoduchší. Máme začiatočnú podmienku. Z nej vypočítame hodnoty v nasledujúcej časovej vrstve. Tieto použijeme na výpočet riešenia v ďalšej časovej vrstve, atď. Na konvergenciu metódy je však potrebné splnenie podmienky na vzťah časového a priestorového kroku. Táto podmienka môže prakticky viesť k nevyhnutnosti zvoliť veľmi malý časový krok - kvôli konvergencii metódy, nie kvôli tomu, že by sme riešenie potrebovali v mnohých tak blízkych časových okamihoch.
V druhom prístupe (implicitná schéma) riešime na každej časovej vrstve sústavu lineárnych rovníc.
E = 50; r = 0.04; D = 0.12; sigma = 0.4;
L = 2;a parametre delenia:
// priestorova premenna n = 20; h = L/n; // casova premenna T = 1; m = 12; k = T/m;
alfa = (r-D)/(sigma^2) - 0.5; beta = (r+D)/2 + (sigma^2)/8 + ((r-D)^2)/(2*sigma^2);
// x = -L, t.j. cena blizka nule function [phi]=phi(tau) phi=0; endfunction // x = L, t.j. cena blizka nekonecnu function [psi]=psi(tau) psi=E*exp(alfa*L + beta*tau).*(exp(L - D*tau) - exp(-r*tau)); endfunctiona začiatočnú podmienku:
// zaciatocna podmienka function [u0]=u0(x) u0=E*exp(alfa*x).*max(0, exp(x)-1); endfunction
ries = zeros(2*n + 1,m+1);a body delenia v čase a v priestore:
x = -L:h:L; tau = 0:k:T;
a = -0.5*(sigma^2)*k/(h^2); // pozdlz diagonaly b = 1 - 2*a; // na diagonalePri výpočte prvej časovej vrstvy (teda hodnoty v čase k) máme nasledovnú pravú stranu:
ps = ries(2:2*n,1)'; ps(1) = ps(1) - a*phi(k); ps(2*n-1) = ps(2*n-1) - a*psi(k);
Môžete použiť funkciu definovanú v súbore gs.sci:
// -------------------------------------------------------------------- // Gauss-Seidelova metoda na riesenie symetrickeho trojdiag.o systemu // -------------------------------------------------------------------- // // gs(a,b,ps,v0,eps) // // 3-diagonalna matica systemu ma: // a - pod a nad diagonalou // b - na diagonale // // ps - prava strana systemu - stlpcovy vektor // v0 - zaciatocna iteracia - stlpcovy vektor // // v - aproximacia riesenia (prepisuje sa v jednotlivych iteraciach) // eps - kriterium pre ukoncenie iteracii: norm(Av-ps)<=eps function [v] = gs(a,b,ps,v0,eps) // matica A pom = length(v0); A(1, 1) = b; for i = 2 : pom A(i, i) = b; A(i, i - 1) = a; A(i - 1, i) = a; end; // cyklus bezi, kym nie je splnena podmienka norm(Av-ps)<=eps v = v0; while norm(A*v - ps) > eps v(1) = (ps(1) - a*v(2))/b; for i = 2:pom-1 v(i) = (ps(i) - a*(v(i-1) + v(i+1)))/b; end v(pom) = (ps(pom) - a*v(pom - 1))/b; end; endfunction
// Navod - pre prvu zlozku: vStare = v(1); vSOR = (ps(1) - a*v(2))/b; v(1) = omega*vSOR + (1-omega)*vStare // omega=1 => povodna SOR // analogicky pre ostatne