How to reduce time in a matlab code?
显示 更早的评论
function [f , t] = jacobi1(n)
tic;
if n < 4
error ('n not in the range')
end
if rem(n,1)~=0
error ('n has to be a whole mumber')
end
n=n^2;
b=zeros(n,1);
b(1,1)=1;
b(n,1)=1;
A = zeros(n,n);
mymatrix=[-1 0 -1 4 -1 0 -1];
for i=1:3
A(i,1:i+3)=mymatrix(5-i:7);
A(n-(3-i),n-(6-i):n)=mymatrix(1:7-i);
end
for i=4:n-3
A(i,i-3:i+3) = mymatrix;
end
epsilon = 1e-3;
f = zeros(n,1);
counter = 0;
flag = 0;
L = tril(A,-1);
D = diag(A);
U = triu(A,1);
B = -1./D.*(L+U);
C = (1./D).*b;
while flag == 0
counter = counter+1;
if counter > 10000
error ('Too many iteration')
end
f_n = (B*f) + C;
if max(abs(f_n-f)/(f_n))<epsilon
flag = 1;
else
f = f_n;
end
end
t=toc;
end
10 个评论
Tzach Berlinsky
2020-12-19
Tzach Berlinsky
2020-12-19
Walter Roberson
2020-12-19
It makes sense for abs(f_n-f)/(f_n) to take a long time, as that is a matrix right division, similar to
abs(f_n-f) * pinv(f_n)
By the way, I suspect you should be taking
max(abs((f_n-f)/f_n))
Tzach Berlinsky
2020-12-19
Walter Roberson
2020-12-19
If you are going to use pinv() then I believe you need
max(abs( (f_n-f) * pinv(f_n) ))
Tzach Berlinsky
2020-12-19
Tzach Berlinsky
2020-12-19
Walter Roberson
2020-12-20
x = B/A solves the system of linear equations x*A = B for x. The matrices A and B must contain the same number of columns.
Mathematically speaking, if A were an invertable matrix and you were using calculations that did not suffer from floating point round-off, then x*A = B could be solved by multiplying both sides on the right by inv(A), giving you x*A*inv(A) = B*inv(A) which is x = B*inv(A) .
And generally speaking, when A is not fully invertable, x = B*pinv(A) can produce a solution -- pinv(A) is, in a way, a generalized inverse that extends to singular matrices as well. For invertable matrices, pinv(A) = inv(A) to within round-off error.
When pinv() is used, then
x_pinv = B*pinv(A)
Breconstructed_pinv = x_pinv * A
sserr_pinv = sum((Breconstructed_pinv - B).^2)
then except for round-off error problems, the sum-squared error between B and the reconstructed B will mathematically be as small as possible.
When / is used:
x_mrd = B/A
Breconstructed_mrd = x_mrd * A
sserr_pinv = sum((Breconstructed_mrd - B).^2)
is not mathematically as small as possible.
However, if I recall correctly, B/A has the property that the reconstructed B will be all 0 except for as many leading entries as there are columns in A -- so in your case with A having one column, the reconstructed B would have one leading non-zero entry and then the rest would be 0.
Using rank() in the situation is just plain WRONG for the situation. It is irrelevant and has nothing to do with the mathematical operation you are trying to perform.
What you are trying to calculate in that step is a measure of how bad the current location is. The calculation (f_n-f)/f_n gives a relative slope, derivatives for each of the components. When the ratio is near 0 in all components, then the location is a good one. If the absolute value of any of the components is too large, then you should keep continuing to look for a better position. But you should not be fitting abs(f_n-f) to f_n as that could give you a large negative component where some element of f_n is negative, and max(abs(f_n-f)/f_n) would pay not attention to the large negative value. Instead, as I have been saying, you should be taking max(abs((f_n-f)/f_n))) as that takes the absolute value of each of the components and looks for the maximum of those. Your code would ignore a component that came out as -1000 if the other components were small enough, whereas my suggestion certainly would not.
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息
2020-12-19
2020-12-20
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
