编写粒子群算法matlab程序，适应度函数为function y=fun(x) y=1/(500^2*x(3)^2*x(1)*x(2)+10/(500^2*40^2*(50*500^2*40^2+1)^2)-(x(3)^2*x(1)*x(2)/500^2*50));

以下是粒子群算法的MATLAB程序，适应度函数为fun：

function [xopt, fopt] = PSO(fun, nvars, lb, ub, MaxIt, nPop, w, c1, c2)
% PSO: Particle Swarm Optimization
%   [xopt, fopt] = PSO(fun, nvars, lb, ub, MaxIt, nPop, w, c1, c2);
%   fun: function handle
%   nvars: number of variables
%   lb: lower bound of variables
%   ub: upper bound of variables
%   MaxIt: maximum number of iterations
%   nPop: population size
%   w: inertial weight
%   c1: cognitive factor
%   c2: social factor
%   xopt: optimal solution
%   fopt: optimal function value

% Initialize the population
pop = repmat(struct('x', [], 'v', [], 'p', [], 'fp', []), nPop, 1);
for i = 1:nPop
    pop(i).x = lb + (ub - lb) .* rand(1, nvars);
    pop(i).v = zeros(1, nvars);
    pop(i).p = pop(i).x;
    pop(i).fp = feval(fun, pop(i).x);
end

% Initialize the best particle and global best
[~, gbest] = min([pop.fp]);
gbest = pop(gbest);
for i = 1:nPop
    if pop(i).fp < gbest.fp
        gbest = pop(i);
    end
end

% Main loop
for t = 1:MaxIt
    % Update the velocity and position of each particle
    for i = 1:nPop
        pop(i).v = w * pop(i).v ...
                  + c1 * rand(1, nvars) .* (pop(i).p - pop(i).x) ...
                  + c2 * rand(1, nvars) .* (gbest.x - pop(i).x);
        pop(i).x = pop(i).x + pop(i).v;
        % Check the boundaries
        pop(i).x = max(pop(i).x, lb);
        pop(i).x = min(pop(i).x, ub);
        % Evaluate the new particle
        pop(i).fp = feval(fun, pop(i).x);
        % Update the personal best
        if pop(i).fp < pop(i).fp
            pop(i).p = pop(i).x;
        end
        % Update the global best
        if pop(i).fp < gbest.fp
            gbest = pop(i);
        end
    end
end

% Output
xopt = gbest.x;
fopt = gbest.fp;
end

使用方法如下：

% Define the fitness function
fun = @(x) 1/(500^2*x(3)^2*x(1)*x(2)+10/(500^2*40^2*(50*500^2*40^2+1)^2)-(x(3)^2*x(1)*x(2)/500^2*50));

% Set the PSO parameters
nvars = 3;          % Number of variables
lb = [0.01 0.01 0.01];  % Lower bounds
ub = [1 1 1];       % Upper bounds
MaxIt = 1000;       % Maximum number of iterations
nPop = 30;          % Population size
w = 1;              % Inertial weight
c1 = 2;             % Cognitive factor
c2 = 2;             % Social factor

% Run the PSO algorithm
[xopt, fopt] = PSO(fun, nvars, lb, ub, MaxIt, nPop, w, c1, c2);

% Display the results
disp('Optimal solution:');
disp(xopt);
disp('Optimal function value:');
disp(fopt);

其中，fun是适应度函数，nvars是变量个数，lb和ub是变量的上下界，MaxIt是最大迭代次数，nPop是种群大小，w、c1和c2是粒子群算法的参数。运行程序即可得到最优解和最优函数值。