I need help with my code

I am having trouble setting up my m-variable for my function c1 = ult_trid(m,b) im supposed to do a forward substitution

but I cannot figure how to declare m.

function Math_5532_6_6_P1

rand('seed',1);
b = rand(4,1);
c = rand(4,1);
d = rand(4,1);
e = rand(4,1);

c(1) = 0;
e(4) = 0;
d1 = d;

[m,d,cri] = lu_fac_trid(c,d,e);

if cri == 0
    disp('A is not invertible');
    return 
end

disp('--------------------------');

disp('max_multiplier');
mm =max(abs(m))
disp('--------------------------');

c1 = ult_trid(m,b);
x = ut_trid(d,e,c1);

%check solution:
r = zeros(4,1);
r(1) = b(1)- (d1(1) * x(1) + e(1) * x(2));
r(2) = b(2)- (c1(2) * x(1) + d1(2) * x(2) + e(2) * x(3));
r(3) = b(3)- (c1(3) * x(2) + d1(3) * x(3) + e(3) * x(4));
r(4) = b(4)- (c1(4) * x(3) + d1(4) * x(4));

disp('The residual vector');
r
disp('--------------------------');

%%%%%%%%%%%%%%%
function [m,d,cri] = lu_fac_trid(c,d,e)
n = length(d);
cri = 1;
m = zeros(n,1);

%Check for invertibility


for k = 1:n-1

    if max(abs(d(k)),abs(c(k+1))) < n*10^(-15)
        cri = 0;

        return
    end

     m(k+1) = c(k+1)/d(k);

     d(k+1) = d(k+1)-m(k+1)*e(k);  
end

if abs(d(n)) < n*10^(-15)
    cri = 0;
    return
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%
function c1 = ult_trid(m,b)
n = length(b);
c1(n) = c1(n)/m(n);
    
    for j=2:n
        c1(j) = c1(j) - m(j)*c1(j);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x=ut_trid(d,e,c1)
    % upper triangle
n=length(d);
x=zeros(n,1);
x(n)=c1(n)/d(n,n);

for i=n-1:-1:1
       
   x(i) = (c1(i)-e(i)*x(i+1)) /d(i);
end
end

end
    

Answer 1

[Walter Roberson]

c1 = ult_trid(m,b);

In your main function, that is the first time you mention c1, and it is an output from ult_trid .

function c1 = ult_trid(m,b)

Okay, within function ult_trid, the inputs are m and b, and the output is c1.

n = length(b);
c1(n) = c1(n)/m(n);

c1 does not exist yet on the right hand side of that assignment. And since we just went through the calling sequence, we know that c1 is not being passed into the function. We also looked at the place that the function was called and have established that c1 does not exist before the function is called, so there is no c1 to pass into the function.

c1(n) = c1(n)/m(n);
    
    for j=2:n
        c1(j) = c1(j) - m(j)*c1(j);
    end

You process c1(n) twice, once explicitly by itself, and the second time when j = n. Meanwhile you never process c1(1) . Are you sure this is wise?

You do not document your functions, so we cannot tell what your intention is in creating and using c1.

As you do not document any "pre-conditions", as outside observers, we must expect that you have accidentally forgotten to check for the possibility that m(n) is 0.

I suspect your code would be a lot easier to write if you used a 2D matrix instead of holding each column as an independent variable.