weixin_47696437 2022-07-11 19:53 采纳率: 85.3%

# matlab求圆心，给出迭代曲线与圆心坐标

###### 问题遇到的现象和发生背景

matlab求圆心，给出迭代曲线与圆心坐标

###### 问题相关代码，请勿粘贴截图
``````%
% clear;
% clc;
% close all;
%
%
% x1 = 556.7058;      y1 = 758.3251;
% x2 = 936.7326;      y2 =1031.4554;
% x3 = 1474.5482;     y3 = 904.8828;
% x4 = 1621.2252;     y4 = 245.3731;
% x5 = 1010.0711;     y5 = 145.4474;
% x6 = 600;     y6 = 427.3114;
%
% x1 = [x1,x2,x3,x4,x5,x6];
% y1 = [y1,y2,y3,y4,y5,y6];
%
% [x0,y0,r] = find_yuanxin(x1,y1);

function [x0,y0,r] = find_yuanxin(x1,y1)
N = 20;
c1 = 2;
c2 = 2;
w = 0.8;
M = 200;
D = 2;

figure(1)
set(gcf,'unit','centimeters','position',[3 12 10.8 10]);
for i = 1:length(x1)-1
line([x1(1,i) x1(1,i+1)],[y1(1,i) y1(1,i+1)]);
hold on
end
line([x1(1,end) x1(1,1)],[y1(1,end) y1(1,1)]);
high = max([max(x1)-min(x1),max(y1)-min(y1)]);
high = ceil(high);
axis([min(x1) min(x1)+high min(y1) min(y1)+high]);

for i = 1:length(x1)
hold on
plot(x1(1,i),y1(1,i),'r.');
end
for i = 1:length(x1)
c = ['D',num2str(i)];
text(x1(i),y1(i),c);
end

xl = min(x1(1),x1(3));    xh = max(x1(1),x1(3));

tt = min([max(x1)-min(x1),max(y1)-min(y1)]);
tx = max(x1) - min(x1);
ty = max(y1) - min(y1);
Scope = [
-1*tx,tx;
-1*ty,ty;
];
Scope1 = [
min(x1)     max(x1);
min(y1)     max(y1);
];

format long;    %设置数据格式
for i = 1:N
x(i,1) = rand*(xh-xl) + xl;   %随机初始化位置
x(i,2) = x(i,1)*(y1(3)-y1(1))/(x1(3)-x1(1)) + (y1(1)*x1(3) - x1(1)*y1(3))/(x1(3)-x1(1));

v(i,1) = rand*(Scope(1,2)-Scope(1,1))+Scope(1,1);   %随机初始化速度
v(i,2) = rand*(Scope(2,2)-Scope(2,1)) + Scope(2,1);
end

for i = 1:N       %利用循环初始化每个粒子的适应度
fp(i) = find_xy(x(i,:),x1,y1);
xp(i,:) = x(i,:);     %每个粒子对应的位置
end
xg = x(1,:);      %给全局最优函数值的位置xg赋初值
for i = 2:N       %通过循环比较找出全局最优函数值的位置xg
if find_xy(x(i,:),x1,y1) < find_xy(xg,x1,y1)
xg = x(i,:);
end
end
for t = 1:M
for i = 1:N
v(i,:) = w*v(i,:)+c1*rand*(xp(i,:)-x(i,:))+c2*rand*(xg-x(i,:));  %更新粒子的速度
x(i,:) = x(i,:) + v(i,:);       %更新粒子的位置

for j = 1:D
if v(i,j) > Scope(j,2);
v(i,j) = Scope(j,2);
elseif v(i,j) < Scope(j,1);
v(i,j) = Scope(j,1);
end
if x(i,j) > Scope1(j,2);
x(i,j) = Scope1(j,2);
elseif x(i,j) < Scope1(j,1);
x(i,j) = Scope1(j,1);
end
end

if find_xy(x(i,:),x1,y1)<fp(i)    %通过比较找出个体最优函数值及其位置
fp(i) = find_xy(x(i,:),x1,y1);
xp(i,:) = x(i,:);
end
if fp(i)<find_xy(xg,x1,y1)        %通过比较找出群体最优函数值及其位置
xg = xp(i,:);
end
end
gbest(t) = find_xy(xg,x1,y1);           %将每个全局最优函数值存储在gbest向量中
end
xmin = xg';      %输出最优位置及函数值
fmin = find_xy(xg,x1,y1);

sita = 0:pi/20:2*pi;
x0 = xmin(1,1);
y0 = xmin(2,1);
r  = fmin*(-1);
plot(x0,y0,'r.');
plot(x0 + r*cos(sita),y0 + r*sin(sita));

end

function r = find_xy(in,X,Y)

for i = 1:length(X)-1
[d(i),dl] = find_dl(in,[X(i),Y(i)],[X(i+1),Y(i+1)]);
end
[d(length(X)),dl] = find_dl(in,[X(end),Y(end)],[X(1),Y(1)]);

r = min(d);
IN = inpolygon(in(1),in(2),[X,X(1)],[Y,Y(1)]);

r = r*(-1);

if IN == 0
r = r + 10000;
end

end

% 求点d1到直线d2-d3的距离，并画线。
function [d,dl] = find_dl(d1,d2,d3)

x1 = d1(1,1);   y1 = d1(1,2);
x2 = d2(1,1);   y2 = d2(1,2);
x3 = d3(1,1);   y3 = d3(1,2);

if x2 == x3
x4 = x1 - 1;
y4 = y1;
elseif y2 == y3;
x4 = x1;
y4 = y1 - 1;
else
k23 = (y3-y2)/(x3-x2);
k1 = -1/k23;
x4 = x1 + 1;
y4 = k1*x4 + y1 - k1*x1;
end

d4 = [x4,y4];
[xx,yy] = node_solve(d1,d4,d2,d3);

dl = [xx,yy];
d = norm(dl-d1);

end

function [X Y] = node_solve( X1,Y1,X2,Y2 )

if X1(1) == Y1(1)
X = X1(1);
k2 = (Y2(2)-X2(2))/(Y2(1)-X2(1));
b2 = X2(2)-k2*X2(1);
Y = k2*X+b2;
end
if X2(1) == Y2(1)
X = X2(1);
k1 = (Y1(2)-X1(2))/(Y1(1)-X1(1));
b1 = X1(2)-k1*X1(1);
Y = k1*X+b1;
end
if X1(1) ~= Y1(1) & X2(1) ~= Y2(1)
k1 = (Y1(2)-X1(2))/(Y1(1)-X1(1));
k2 = (Y2(2)-X2(2))/(Y2(1)-X2(1));
b1 = X1(2)-k1*X1(1);
b2 = X2(2)-k2*X2(1);
if k1 == k2
X = [];
Y = [];
else
X = (b2-b1)/(k1-k2);
Y = k1*X+b1;
end
end

end

``````

###### 我想要达到的结果

matlab求圆心，给出迭代曲线与圆心坐标

• 写回答

#### 1条回答默认 最新

• Wayne_Fine Matlab领域优质创作者 2022-07-11 22:57
关注

圆心坐标直接在窗口里就能查看啊（中间那个红点）：

加上命令行输出以后：

``````%
clear;
clc;
close all;

x1 = 556.7058;      y1 = 758.3251;
x2 = 936.7326;      y2 =1031.4554;
x3 = 1474.5482;     y3 = 904.8828;
x4 = 1621.2252;     y4 = 245.3731;
x5 = 1010.0711;     y5 = 145.4474;
x6 = 600;     y6 = 427.3114;

x1 = [x1,x2,x3,x4,x5,x6];
y1 = [y1,y2,y3,y4,y5,y6];

[x0,y0,r] = find_yuanxin(x1,y1);

function [x0,y0,r] = find_yuanxin(x1,y1)
N = 20;
c1 = 2;
c2 = 2;
w = 0.8;
M = 200;
D = 2;

figure(1)
set(gcf,'unit','centimeters','position',[3 12 10.8 10]);
for i = 1:length(x1)-1
line([x1(1,i) x1(1,i+1)],[y1(1,i) y1(1,i+1)]);
hold on
end
line([x1(1,end) x1(1,1)],[y1(1,end) y1(1,1)]);
high = max([max(x1)-min(x1),max(y1)-min(y1)]);
high = ceil(high);
axis([min(x1) min(x1)+high min(y1) min(y1)+high]);

for i = 1:length(x1)
hold on
plot(x1(1,i),y1(1,i),'r.');
end
for i = 1:length(x1)
c = ['D',num2str(i)];
text(x1(i),y1(i),c);
end

xl = min(x1(1),x1(3));    xh = max(x1(1),x1(3));

tt = min([max(x1)-min(x1),max(y1)-min(y1)]);
tx = max(x1) - min(x1);
ty = max(y1) - min(y1);
Scope = [
-1*tx,tx;
-1*ty,ty;
];
Scope1 = [
min(x1)     max(x1);
min(y1)     max(y1);
];

format long;    %设置数据格式
for i = 1:N
x(i,1) = rand*(xh-xl) + xl;   %随机初始化位置
x(i,2) = x(i,1)*(y1(3)-y1(1))/(x1(3)-x1(1)) + (y1(1)*x1(3) - x1(1)*y1(3))/(x1(3)-x1(1));

v(i,1) = rand*(Scope(1,2)-Scope(1,1))+Scope(1,1);   %随机初始化速度
v(i,2) = rand*(Scope(2,2)-Scope(2,1)) + Scope(2,1);
end

for i = 1:N       %利用循环初始化每个粒子的适应度
fp(i) = find_xy(x(i,:),x1,y1);
xp(i,:) = x(i,:);     %每个粒子对应的位置
end
xg = x(1,:);      %给全局最优函数值的位置xg赋初值
for i = 2:N       %通过循环比较找出全局最优函数值的位置xg
if find_xy(x(i,:),x1,y1) < find_xy(xg,x1,y1)
xg = x(i,:);
end
end
for t = 1:M
for i = 1:N
v(i,:) = w*v(i,:)+c1*rand*(xp(i,:)-x(i,:))+c2*rand*(xg-x(i,:));  %更新粒子的速度
x(i,:) = x(i,:) + v(i,:);       %更新粒子的位置

for j = 1:D
if v(i,j) > Scope(j,2);
v(i,j) = Scope(j,2);
elseif v(i,j) < Scope(j,1);
v(i,j) = Scope(j,1);
end
if x(i,j) > Scope1(j,2);
x(i,j) = Scope1(j,2);
elseif x(i,j) < Scope1(j,1);
x(i,j) = Scope1(j,1);
end
end

if find_xy(x(i,:),x1,y1)<fp(i)    %通过比较找出个体最优函数值及其位置
fp(i) = find_xy(x(i,:),x1,y1);
xp(i,:) = x(i,:);
end
if fp(i)<find_xy(xg,x1,y1)        %通过比较找出群体最优函数值及其位置
xg = xp(i,:);
end
end
gbest(t) = find_xy(xg,x1,y1);           %将每个全局最优函数值存储在gbest向量中
end
xmin = xg';      %输出最优位置及函数值
fmin = find_xy(xg,x1,y1);

sita = 0:pi/20:2*pi;
x0 = xmin(1,1)   %圆心横坐标
y0 = xmin(2,1)   %圆心纵坐标
r  = fmin*(-1);
plot(x0,y0,'r.');
plot(x0 + r*cos(sita),y0 + r*sin(sita));

end

function r = find_xy(in,X,Y)

for i = 1:length(X)-1
[d(i),dl] = find_dl(in,[X(i),Y(i)],[X(i+1),Y(i+1)]);
end
[d(length(X)),dl] = find_dl(in,[X(end),Y(end)],[X(1),Y(1)]);

r = min(d);
IN = inpolygon(in(1),in(2),[X,X(1)],[Y,Y(1)]);

r = r*(-1);

if IN == 0
r = r + 10000;
end

end

% 求点d1到直线d2-d3的距离，并画线。
function [d,dl] = find_dl(d1,d2,d3)

x1 = d1(1,1);   y1 = d1(1,2);
x2 = d2(1,1);   y2 = d2(1,2);
x3 = d3(1,1);   y3 = d3(1,2);

if x2 == x3
x4 = x1 - 1;
y4 = y1;
elseif y2 == y3;
x4 = x1;
y4 = y1 - 1;
else
k23 = (y3-y2)/(x3-x2);
k1 = -1/k23;
x4 = x1 + 1;
y4 = k1*x4 + y1 - k1*x1;
end

d4 = [x4,y4];
[xx,yy] = node_solve(d1,d4,d2,d3);

dl = [xx,yy];
d = norm(dl-d1);

end

function [X Y] = node_solve( X1,Y1,X2,Y2 )

if X1(1) == Y1(1)
X = X1(1);
k2 = (Y2(2)-X2(2))/(Y2(1)-X2(1));
b2 = X2(2)-k2*X2(1);
Y = k2*X+b2;
end
if X2(1) == Y2(1)
X = X2(1);
k1 = (Y1(2)-X1(2))/(Y1(1)-X1(1));
b1 = X1(2)-k1*X1(1);
Y = k1*X+b1;
end
if X1(1) ~= Y1(1) & X2(1) ~= Y2(1)
k1 = (Y1(2)-X1(2))/(Y1(1)-X1(1));
k2 = (Y2(2)-X2(2))/(Y2(1)-X2(1));
b1 = X1(2)-k1*X1(1);
b2 = X2(2)-k2*X2(1);
if k1 == k2
X = [];
Y = [];
else
X = (b2-b1)/(k1-k2);
Y = k1*X+b1;
end
end

end
``````

打印结果：
x0 =
1.105632284006582e+03
y0 =
5.736136195988425e+02

本回答被题主选为最佳回答 , 对您是否有帮助呢?
评论 编辑记录

• 系统已结题 7月26日
• 已采纳回答 7月18日
• 创建了问题 7月11日

#### 悬赏问题

• ¥15 python结合Matlab仿真忆阻器
• ¥35 有人会注册whatsaop协议号吗？