# I need help with my code

2 views (last 30 days)
daniel choudhry on 20 Oct 2020
Answered: Walter Roberson on 20 Oct 2020
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

Walter Roberson on 20 Oct 2020
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.