I need help with my code

Question

daniel choudhry 2020-10-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/620543-i-need-help-with-my-code

回答： Walter Roberson 2020-10-20

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
    

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2020-10-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/620543-i-need-help-with-my-code#answer_519633

在 MATLAB Online 中打开

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.