Why do I keep getting a zero for x(9) for Backward subsitution?

This is my code for backward subsitution. I keep getting zero for x(9) each run. I'm not sure what I'm doing wrong.

%% b
 U = triu(randi([0,20],10,10));
 b = randi([0,20],10,1);
 x = backward_sub(U,b)

x = 10×1

0.1875 2.7143 1.4444 3.1667 1.0000 1.3333 Inf 2.2857 0 0.5000

 %%
 function [x] = backward_sub(U,b) %Initializing function
    n= length(b); % Gets the length of b
    x = zeros(n,1);
    x(n) = b(n) / U(n,n); % gets the known value 
    for i = n-1:-1:1 % loops backward over the compents of U
        sum = 0; %intializing the sum
        for j = 1:i-1 %loops over the coments of the sum
            sum = sum + (U(i,j)*x(j)); % takes the sum of the components
        end
        x(j) = (b(j) - sum) / U(i,i); % finds x(j) for every fixed k. 
    end
 end

Answer 1

[Torsten]

Don't generate random numbers that can become 0 on the diagonal. Or use your code to check if all U(i,i) are not equal 0.

Further see the changes made in backward_sub.

%% b
 U = triu(randi([1,20],10,10))

U = 10×10

17 10 15 15 4 11 17 6 3 8 0 2 14 17 20 6 2 16 8 17 0 0 20 17 5 19 11 11 19 19 0 0 0 8 8 5 17 15 14 14 0 0 0 0 12 10 19 2 7 11 0 0 0 0 0 5 14 4 6 9 0 0 0 0 0 0 5 1 20 7 0 0 0 0 0 0 0 11 6 15 0 0 0 0 0 0 0 0 6 7 0 0 0 0 0 0 0 0 0 17

 b = randi([0,20],10,1);
 x = backward_sub(U,b)

x = 10×1

-21.8157 46.7439 -12.5169 3.8361 -2.1406 11.4292 -4.2927 -0.4385 1.4216 0.3529

 U*x-b

ans = 10×1

1.0e-13 * 0.6040 -0.1066 -0.5418 0 0.0355 0.0622 -0.0178 0 0 0

 %%
 function [x] = backward_sub(U,b) %Initializing function
    n= length(b); % Gets the length of b
    x = zeros(n,1);
    x(n) = b(n) / U(n,n); % gets the known value 
    for i = n-1:-1:1 % loops backward over the compents of U
        sum = 0; %intializing the sum
        for j = i+1:n %loops over the coments of the sum
            sum = sum + (U(i,j)*x(j)); % takes the sum of the components
        end
        x(i) = (b(i) - sum) / U(i,i); % finds x(j) for every fixed k. 
    end
 end

[Jessica Arroyo]

I have to have U in a range between 0 and 20.

[Torsten]

U*x-b what is this part for?

To check if your solution is correct. It should come out to be the zero vector, shouldn't it ?

I have to have U in a range between 0 and 20.

As I said: If you must have U in the range 0 to 20, it can happen that some elements on the diagonal of U become zero. Then your matrix U is rank deficient, and usual back-substitution doesn't work.

[Jessica Arroyo]

Why did the for loop for j change. do you mind explaining?

[Torsten]

You must take the terms to the right of the diagonal to the right-hand side of the i-th equation and then divide by the diagonal element U(i,i) to get x(i). That's why you first sum the elements to the right of the diagonal (from i+1 to n), subtract them from the right-hand side b(i) (b(i)-sum) and then divide by the diagonal element U(i,i).

You should just right down the programmed operations - then you'll see what's going on.

[Jessica Arroyo]

Ohh I see, thank you so much! :-)

Answer 2

[John D'Errico]

First: NEVER use the variable name sum as a variable name. Why not? Becaause one time you will want to use the FUNCTION sum. And then you will post an anquished question on Answers asking why sum no longer works for you.

As far as why your code does not work? It has a bug in it. Correctly written code will work. ;-)

U = triu(randi([0,20],10,10));
b = randi([0,20],10,1);

What solution do we expect?

U\b

ans = 10×1

-0.9462 2.0157 1.3474 -5.4897 -0.1753 1.2635 0.4722 -0.2500 -0.4286 1.0000

Does my corrected version work?

x = backward_sub(U,b)

x = 10×1

-0.9462 2.0157 1.3474 -5.4897 -0.1753 1.2635 0.4722 -0.2500 -0.4286 1.0000

You might see there are quite a few changes I made. Two spelling changes in comment, because you were not spelling the word component correctly, and it bugged me. I also refused to name a variable sum. But more importantly, I used the correct loops. And you seem to have been using j in a few places where you needed to use i.

function [x] = backward_sub(U,b) %Initializing function
    n= length(b); % Gets the length of b
    x = zeros(n,1);
    x(n) = b(n) / U(n,n); % gets the known value 
    for i = n-1:-1:1 % loops backward over the components of U
        sumi = 0; %intializing the sum
        for j = i+1:n %loops over the components of the sum
            sumi = sumi + (U(i,j)*x(j)); % takes the sum of the components
        end
        x(i) = (b(i) - sumi) / U(i,i); % finds x(i)
    end
 end

[Jessica Arroyo]

Thank you!