How to vectorize this code?
显示 更早的评论
Can anyone please help me to write a vectorized version of following code assignment. I tried by using meshgrid instead of two outer loops of x and y values:
=============================
clc
clear all
close all
r1=[0;0];
r2=[1;1];
r3=[-1;-1];
for x=-2:0.03:2
for y=-2:0.03:2
X=[x;y];
for j=1:30
f=[X(1)^3 - X(2), X(2)^3 - X(1)]';
Jf=[3*X(1)^2, -1; -1, 3*X(2)^2];
X=X-Jf\f;
end
format long
X;
if norm(X-r1)<1e-8 % if the distance between X and r1 is less than 10^-8
c='b'; % use the color red
elseif norm(X-r2)<1e-8
c='m'; % use the color blue
elseif norm(X-r3)<1e-8
c='r'; % use the color green
else % if not close to any of the roots
c='y'; % use the color black
end
plot(x,y,'.','Color',c,'Markersize',40);
hold on
end
end
=============================================
Thanx
采纳的回答
Roger Stafford
2013-9-1
Your outer two for-loops can be vectorized, but the innermost one which involves the iteration process cannot - at least as far as I am aware.
[X,Y] = meshgrid(-2:0.03:2);
for j = 1:30
D = 9*X.^2.*Y.^2-1; % Avoid matrix division
N = 6*X.^2.*Y.^3+2*X.^3; % Use direct expression for newton-raphson formula
X = (6*X.^3.*Y.^2+2*Y.^3)./D;
Y = N./D;
tr = (X-0).^2+(Y-0).^2<1e-16; % Color the dots red for true tr
tb = (X-1).^2+(Y-1).^2<1e-16; % Color the dots blue for true tb
tg = (X+1).^2+(Y+1).^2<1e-16; % Color the dots green for true tg
ty = ~(tr|tb|tg); % Color the other dots black
end
plot(X(tr),Y(tr),'r.',etc....
By the way, there are six more roots to your equations. They are all complex numbers on the unit circle in the complex plane. For example, (x,y) = (1i,-1i) is a root.
3 个评论
Athar
2013-9-1
Roger Stafford
2013-9-1
Don't try to plot the black points. They are far, far out because the newton-raphson algorithm is not converging for them but rather diverging. That is true for both the vectorized method and your nested for-loop method. It only works when the initial points are within some reasonable neighborhood of the solution. To attempt a plot of all the diverging points will cause the scaling of the plot to be so huge as to make all valid solution points look as if they are packed right at the origin. Look at the values in your plot axes to see this. The last few lines should be:
tr = (X-0).^2+(Y-0).^2<1e-16; % Color the dots red for true tr
tb = (X-1).^2+(Y-1).^2<1e-16; % Color the dots blue for true tb
tg = (X+1).^2+(Y+1).^2<1e-16; % Color the dots green for true tg
end
plot(X(tr),Y(tr),'ro',X(tb),Y(tb),'bo',X(tg),Y(tg),'go')
Roger Stafford
2013-9-1
I would think it would be far more informative to plot the original positions, rather than final positions, of the points which converged to one of your three solutions. You already know where the final points are to be found but you don't know where they came from.
You can do this by saving a copy of the original X and Y and using the 't' logical arrays with them. This will make plots of the "neighborhoods" I mentioned previously.
[X,Y] = meshgrid(-2:0.03:2);
X0 = X; Y0 = Y;
.....
plot(X0(tr),Y0(tr),'r*',X0(tb),Y0(tb),'b*',X0(tg),Y0(tg),'g*')
更多回答（1 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
