Forward and backward substitution of triangular matricies

I am trying to code a forward and backward substitution for a homework problem.

Here is my forward code with an example.

function x=myForwardSubstitution(L,b)

d=size(L,1);

L=[1,0;1,1];

b=[1;1];

for i=1:d

for j=1:i-1

b(i)=b(i)-L(i,j)*b(j);

end

b(i)=b(i)/L(i,i);

end

x=b;

end

and here is the backward one also with an example.

function x=myBackwardSubstitution(U,b)

d=size(U,1);

U=[1,2;0,1];

b=[1;1];

for i=d:-1:1

for j=1:i-1

b(i)=b(i)-U(i,j)*b(j);

end

b(i)=b(i)/U(i,i);

end

x=b;

end

I have underlined the lines which get an error. The error says "Index exceeds the number of array elements (2)" Have I defined my i and j wrong?

Thanks in advance

Answer 1

[Abdullah Zaihi]

clear all

clc

n = 50; % change it to 5000

% generate a set of n orthogonal vectors to use as eigenvectors

B = rand(n, n);

[Q, ~] = qr(B); %QR decomposition Q(n,n) R(n,m) A(n,n)=QR the R is niglicted

% generate a set of n eigenvalues between 1 <= lambda <= 1000

lambda = linspace(1, 1000, n);

D = diag(lambda);

% generate the SPD dense matrix explicitly

A = Q * D * Q';

% compute the Cholesky decomposition

R = chol(A); % matlab returns the *upper* triangular factor A=R'*R A suppose to be symmetric matrix mxm

% YOUR CODE GOES HERE

T=R'*R; %verification T=A

b=rand(n,1);

RT=R';

o=b; %for not messing with b vector

%forward substituation

y(1)=o(1)/RT(1,1);

for i=1:n

for j=1:i-1

o(i)=o(i)-RT(i,j)*y(j);

end

y(i,1)=o(i)/RT(i,i);

end

z=RT\b; %verify the forward substituation z=y

w=y; %for not messing with y vector

%back substituation

x(n,1) = w(n) / R(n,n);

for i=n-1:-1:1

for j=i+1:n

w(i)=w(i)-R(i,j)*x(j);

end

x(i) = w(i)/R(i,i);

end

v=R\y; %verify the back substituation v=x

r=norm(A*x-b) %checking the residual for forward and back sub.

m=A\b; %verify using backslash method m=v=x

r_m=norm(A*m-b) %checking the residual for backslash method

Answer 2

[Steven Lord]

function x=myForwardSubstitution(L,b)
    d=size(L,1);
   L=[1,0;1,1];

You define d using the L matrix with which the user of your code called your function, but then you throw away the user's input and replace it with a matrix that may have fewer rows than that input. As a similar but simpler example, it's as though you ran these four lines of code.

L = [1, 2, 3; 4, 5, 6; 7, 8, 9];
d = size(L, 1)
L = [1, 0; 1, 1];
L(d, 1) % will error as L no longer has 3 rows

I would eliminate the code in your function where you overwrite the user's input, or move that code outside the function and pass the result of running that code into your function.