Can someone help me figure out what is wrong with my function c1 = ult_trid(m,b) ??

Question

daniel choudhry 2020-10-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/608761-can-someone-help-me-figure-out-what-is-wrong-with-my-function-c1-ult_trid-m-b

回答： James Tursa 2020-10-9

I cannot get my code to run because it tells me:Index in position 2 exceeds array bounds (must not exceed 1).

which i cannot figure out what part of my function c1 = ult_trid(m,b) is wrong. can someone help me?

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) = d(k);
        d(k+1) = c(k+1);
end
if abs(d(n)) < n*10^(-15)
    cri = 0;
    return
end
 for i=n:-1:k+1
        factor= c(i)/c(k);
        factor2= e(i)*d(k);
       m(i)=c(i)-c(k)*(factor);
       
        d(i)=d(i)-d(k)*(factor2);
 end
end
%%%%%%%%%%%%%%%%%%%%%%%%%
function c1 = ult_trid(m,b)
    
    [brows, bcol] = size(b);
[mrows, mcol] = size(m);
c1 = b;
for k = 1:mrows
    for j = 1:k-1
        c1(k) = c1(j) - m(j)*c1(j);
    end
    c1(k) = c1(k)/m(k);
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