LU dekompozicija trodijagonalne matrice

LUtri[n_, a_, b_, c_, y_] := 
      Module[{d, e, f, z, x}, 
 (* 
          Treba unijeti red matrice n,  vektore glavne a,  
          gornje b i donje c dijagonale,  
          te vektor slobodnih koeficijenata y *) 
 
        d =   e =  f = z =  x = 
                  Table[0, {i, 
                      n}]; 
                                            (* 
          Izracunavanje dekompozicije *) 
e[[1]] = 
          a[[1]]; f[[1]] = b[[1]]; 
        d[[1]] = c[[1]]/e[[1]];  
        Do[e[[i]] = a[[i]] - d[[i - 1]] f[[i - 1]]; 
          f[[i]] = b[[i]];  
          d[[i]] = c[[i]]/e[[i]], {i, 2, 
            n - 1}]; 
e[[n]] = 
          a[[n]] - 
            d[[n - 1]] f[[n - 
                    1]]; 
                                            (* 
          Rjesavanje sustava LUx = y *) 
z[[1]] = 
          y[[1]]; 
Do[
          z[[i]] = y[[i]] - d[[i - 1]] z[[i - 1]], {i, 2, 
            n}]; 
x[[n]] = 
          z[[n]]/e[[n]]; 
 
        Do[x[[n - i]] = (z[[n - i]] - 
                  f[[n - i]] x[[n - i + 1]])/e[[n - i]], {i, 
            n - 1}]; 
x // N]

Definiranje trodijagonalne matrice i poziv modula LUtri

n = 5; 
a = Table[2  i, {i, n}];  b = Table[i - 1, {i, n - 1}];  
    c = Table[i + 1, {i, n - 1}];  
    y = Table[i, {i, n}];
A = Table[0, {i, n}, {j, n}]; 
Do[A[[i, i]] = a[[i]], {i, n}];
Do[A[[i, i + 1]] = b[[i]]; 
        A[[i + 1, i]] = c[[i]], {i, 
          n - 1}];
Print["A=", MatrixForm[A],  "   y=", 
      MatrixForm[y]]
x = LUtri[n, a, b, c, y];
Print["x = ", x ]

Izvod formula

n = 5;
a = {a1, a2, a3, a4, a5};    b = {b1, b2, b3, 
        b4};            c = {c1, c2, c3, 
        c4};
d = {d1, d2, d3, d4};            e = {e1, e2, e3, e4, 
        e5};     f = {f1, f2, f3, f4};
A = L = U = 
            Table[0, {i, n}, {j, n}]; 
Do[A[[i, i]] = a[[i]]; L[[i, i]] = 1; 
        U[[i, i]] = e[[i]], {i, n}];
Do[A[[i, i + 1]] = b[[i]]; A[[i + 1, i]] = c[[i]];  
        L[[i + 1, i]] = d[[i]]; 
        U[[i, i + 1]] = f[[i]], {i, 
          n - 1}];
Print["A=", MatrixForm[A], "   L=", 
      MatrixForm[L], "    U=", MatrixForm[U]]
Print["  LU=", MatrixForm[L . U] ]

Programiranje

\(\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ 
      Podaci\ *) \)\(
\)\(\(n = 5;\)

    a = Table[2  i, {i, n}]; \ \ \ \ \ \ b = 
      Table[i - 1, {i, n - 1}]; \ \ \ \ c = 
      Table[i + 1, {i, n - 1}];

    \(y = Table[i, {i, n}];\)

    \(mat = Table[0, {i, n}, {j, n}];\)

    \(Do[mat[\([i, i]\)] = a[\([i]\)], {i, n}];\)

    Do[mat[\([i, i + 1]\)] = b[\([i]\)]; 
      mat[\([i + 1, i]\)] = c[\([i]\)], {i, n - 1}]

    Print["A=", MatrixForm[mat], "   y=", 
      MatrixForm[y]]\)\)\)

\(\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ 
      LU - dekompozicija\ *) \)\(
\)\(\(d = \(e = \(f = 
            Table[0, {i, n}]\)\);\)

    e[\([1]\)] = a[\([1]\)]; \ f[\([1]\)] = b[\([1]\)]; \ 
    d[\([1]\)] = c[\([1]\)]/e[\([1]\)];

    \(Do[e[\([i]\)] = 
          a[\([i]\)] - 
            d[\([i - 1]\)] f[\([i - 1]\)]; 
f[\([i]\)] = 
          b[\([i]\)]; 
d[\([i]\)] = 
          c[\([i]\)]/e[\([i]\)], 
{i, 2, 
          n - 1}];\)

    \(\(e[\([n]\)] = 
        a[\([n]\)] - d[\([n - 1]\)] f[\([n - 1]\)];\)\(
\)
    \)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ Ispis\ \ *) 

    \(L = \(U = Table[0, {i, n}, {j, n}]\);\)

    \(Do[L[\([i, i]\)] = 1; 
        U[\([i, i]\)] = e[\([i]\)], {i, n}];\)

    \(Do[\ L[\([i + 1, i]\)] = d[\([i]\)]; 
        U[\([i, i + 1]\)] = f[\([i]\)], {i, n - 1}];\)

    \(\(Print["L=", MatrixForm[L], "  U=", 
      MatrixForm[U]]\)\(
\)
    \)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ Kontrola\ *) 

    Print["A=", MatrixForm[L . U]]\)\)\)

z = x = 
          Table[0, {i, n}];

    z[[1]] = y[[1]];

    Do[z[[i]] = y[[i]] - d[[i - 1]] z[[i - 1]], {i, 2, 
          n}];

    Print["z=", z // N]

    x[[n]] = z[[n]]/e[[n]];

     Do[
        x[[n - i]] = (z[[n - i]] - 
                f[[n - i]] x[[n - i + 1]])/e[[n - i]], {i, 
          n - 1}];

    Print["x=", x // N]

Test primjer

n = 3000;  SeedRandom[13];
a = Table[Random[Real, {0, 100}], {i, n}];            b = 
      Table[Random[Real, {0, 100}], {i, n - 1}];
c = Table[Random[Real, {0, 100}], {i, n - 1}]; 
    y = Table[Random[Real, {0, 100}], {i, n}];

\(\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ 
      LU\ dekompozicija\ *) \)\(
\)\(t1 = \(Timing[
x = LUtri[n, a, b, c, y]; 
 (*\ 
            Print["   x = ", 
              x // N]\ *) 
]\)[\([1]\)]\)\)\)

\(\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ 
      Mathematica\ *) \)\(
\)\(\ \)\(\(A = 
        Table[0, {i, n}, {j, n}];\)\ 

    \(Do[A[\([i, i]\)] = a[\([i]\)], {i, n}];\)

    \(\(Do[A[\([i, i + 1]\)] = b[\([i]\)]; 
        A[\([i + 1, i]\)] = c[\([i]\)], {i, n - 1}];\)\(
\)
    \) (*\ Print["A=", MatrixForm[A], \ "   y=", 
        MatrixForm[y]]\ \ *) 

    \(\(x2 = Timing[LUBackSubstitution[LUDecomposition[A], y]]\) \); 
    x2[\([1]\)]\)\(\ \)\)\)