I need help with my code

1 次查看（过去 30 天）
daniel choudhry
daniel choudhry 2020-10-20
回答： 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

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

采纳的回答

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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by