force-direct力导向算法 matlab实现

简单实现力导向算法，节点不要太多。
设置参数需根据节点个数情况而定，多次尝试效果较好

clc;clear;close all
%% 导入
% connect_list=[1 2;2 3;3 4;4 5;5 1]
load('connect_list.mat')
%% 参数设定
kr=10;
ks=0.0001;
x0=300;
dt=1;
maxmove=50;
dd=200;
%% 节点计算
point=[connect_list(:,1);connect_list(:,2)];
point=sort(point);
point_save(1,1)=point(1,1);
n=1;
for i=2:length(point)
    if point(i)~=point_save(n)
        n=n+1;
        point_save(n,1)=point(i,1);
    end
end
point_connect=point_save;
for i=1:length(point_connect(:,1))
    
    for j=1:length(connect_list(:,1))
        if connect_list(j,1)==point_connect(i,1)
            n=2;
            for k=2:length(point_connect(1,:))
                if point_connect(i,k)~=0
                    n=n+1;
                end
            end
            point_connect(i,n)=connect_list(j,2);
        else if connect_list(j,2)==point_connect(i,1)
                n=2;
                for k=2:length(point_connect(1,:))
                    if point_connect(i,k)~=0
                        n=n+1;
                    end
                end
                point_connect(i,n)=connect_list(j,1);
            end
        end
    end
end

for i=1:length(point_save)
    point_save(i,2)=rand(1)*100;
    point_save(i,3)=rand(1)*100;
    point_save(i,4)=rand(1)*100;
end
for n=1:dd
%% 洛伦兹
force1=zeros(length(point_save(:,1)),4);
for i=1:length(point_save(:,1))
    for j=i+1:length(point_save(:,1))
        dxyz=point_save(i,2:4)-point_save(j,2:4);
        dsq=dxyz.^2;
        d=sqrt(dsq(1)+dsq(2)+dsq(3));
        fxyz=kr/d^2*dxyz./d;
        force1(i,2:4)=force1(i,2:4)-fxyz;
        force1(j,2:4)=force1(j,2:4)+fxyz;
    end
end

%% 弹簧力
force2=zeros(length(point_save(:,1)),4);
for i=1:length(point_connect(:,1))
     for j=i+1:length(point_save(:,1))
        dxyz=point_save(i,2:4)-point_save(j,2:4);
        dsq=dxyz.^2;
        d=sqrt(dsq(1)+dsq(2));
        fxyz=ks*(d-x0)*dxyz./d;
        force2(i,2:4)=force2(i,2:4)+fxyz;
        force2(j,2:4)=force2(j,2:4)-fxyz;
    end
end
%% 力转化位移距离
force=force1+force2;
distence=force.*dt;
% distence(;,1)=point_save(:,1);
for i=1:length(distence(:,1))
    if distence(i,2)>maxmove
        distence(i,2)=maxmove;
    else if distence(i,2)<-maxmove
            distence(i,2)=-maxmove;
        end
    end
    if distence(i,3)>maxmove
        distence(i,3)=maxmove;
    else if distence(i,3)<-maxmove
            distence(i,3)=-maxmove;
        end
    end
        if distence(i,4)>maxmove
        distence(i,4)=maxmove;
    else if distence(i,4)<-maxmove
            distence(i,4)=-maxmove;
        end
    end
end
point_save=point_save+distence;
n
end
%% 画图
for i=1:length(connect_list(:,1))
    plot3([point_save(connect_list(i,1),2),point_save(connect_list(i,2),2)],...
        [point_save(connect_list(i,1),3),point_save(connect_list(i,2),3)],...
        [point_save(connect_list(i,1),4),point_save(connect_list(i,2),4)],...
        '-','color',[142,154,175]/255);hold on
    plot3(point_save(connect_list(i,1),2),point_save(connect_list(i,1),3),point_save(connect_list(i,1),4),...
        'ro','LineWidth',1,'MarkerSize', 2,'MarkerEdgeColor',[87,96,105]/255,...
        'markerfacecolor', [247,68,97]/255)
    hold on
end
% 如果点的标号是无序的用这个画图
% for i=1:length(connect_list(:,1))
%     for j=1:length(point_save(:,1))
%         if point_save(j,1)==connect_list(i,1)
%             p1=j;break
%         end
%     end
%     for j=1:length(point_save)
%         if point_save(j,1)==connect_list(i,2)
%             p2=j;break
%         end
%     end
%     plot3([point_save(p1,2),point_save(p2,2)],...
%         [point_save(p1,3),point_save(p2,3)],...
%         [point_save(p1,4),point_save(p2,4)],...
%         '-','color',[142,154,175]/255);hold on
%     plot3(point_save(p1,2),point_save(p1,3),point_save(p1,4),...
%         'ro','LineWidth',1,'MarkerSize', 2,'MarkerEdgeColor',[87,96,105]/255,...
%         'markerfacecolor', [247,68,97]/255)
%     hold on
% end