How to minimise potential energy using MATLAB?

Question

Given m ∈ N of immobile electric charges at some points of the plane and n ∈ N number of further charged particles along a circle centered at the origin, in such a way that the sum potential of the system is minimized. Charges are equal. How to place n charges so potential would be minimal? Which algorithm should be used? I have this code which calculates total potential energy:

function example = potential_energy()

    clf;

    % given charged points
    % all charges are equal

    % m number of immobile charges in random positions
    m = 5;
    % n number of charges along the circle with fixed radius
    n = 10;

    % task is minimize total potential energy

    % generate and plot immobile charges
    mx=rand(1,m)*10;
    my=rand(1,m)*10;
    plot(mx,my,"b*");
    hold on;

    % generate and plot charges along the circle
    angles = linspace(0, 2*pi, 720);
    radius = 5;
    xCenter = 5;
    yCenter = 5;
    x = radius * cos(angles) + xCenter; 
    y = radius * sin(angles) + yCenter;
    randomIndexes = randperm(length(angles), n);
    nx = x(randomIndexes);
    ny = y(randomIndexes);
    plot(nx, ny, 'r*');

    % join m and n charges into one array
    x = [mx,nx];
    y = [my,ny];

    total_potential = 0;
    % using this loop, total potential is calculated
    for i = 1:length(x)
        for j = i+1:length(x)
            distance = sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
            % In order to calculate potential energy: U=kQ/r
            % All charges are equal so for simplicity let's take Q(charge) as 1/k
            % So potential energy can be calculated using U=1/r 
            potential = 1/distance;
            total_potential = total_potential + potential;
        end
    end

    disp(total_potential);

end

How total potential energy can be minimized? Which algorithm should be used?

Answer 1

You can use standard tools from the MATLAB optimization toolbox like fminsearch or fminunc.

Just put your function to evaluate the potential in a MATLAB file, e.g. test_potE.m

function total_potential = test_potE(angles,mx,my,radius,xCenter,yCenter)
% angles ... variables to be optimized
% the rest of the input variables are fixed but required to compute the output total_potential 

 nx = radius * cos(angles) + xCenter; 
 ny = radius * sin(angles) + yCenter;

 % join m and n charges into one array
 x = [mx,nx];
 y = [my,ny];

 total_potential = 0;
 % using this loop, total potential is calculated
 for i = 1:length(x)
    for j = i+1:length(x)
        distance = sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
        % In order to calculate potential energy: U=kQ/r
        % All charges are equal so for simplicity let's take Q(charge) as 1/k
        % So potential energy can be calculated using U=1/r 
        potential = 1/distance;
        total_potential = total_potential + potential;
    end
 end
 ```

Then, optimize it in the main file. For your problem it doesn't matter if constraint or unconstraint optimization tools are used. In case local minima exist you have to repeat the optimization several times to find the global one.

clf;

% given charged points
% all charges are equal

% m number of immobile charges in random positions
m = 5;
% n number of charges along the circle with fixed radius
n = 10;

% task is minimize total potential energy

% generate and plot immobile charges
mx=rand(1,m)*10;
my=rand(1,m)*10;
plot(mx,my,'b*');
hold on;

% generate and plot charges along the circle
radius = 5;
xCenter = 5;
yCenter = 5;

% the variables to be optimized are the angles on the circle
% initialized as given above
angles = linspace(0, 2*pi, 720);
randomIndexes = randperm(length(angles), n);
angles = angles(randomIndexes);

% plot the initial points on the circle
nx = radius * cos(angles) + xCenter; 
ny = radius * sin(angles) + yCenter;
plot(nx, ny, 'r*');

% we need to specify the parameters of the function that 
% computes the potential, because fminsearch,... only takes the function handle
% and the initial angles
f = @(angles)test_potE(angles,mx,my,radius,xCenter,yCenter);

total_potential = f(angles);
disp(total_potential);

% .. display stuff during learning ...
options = optimset('Display','iter','PlotFcns',@optimplotfval);

% optimization
angles1 = fminsearch(f,angles,options) 
% ... or ...
angles2 = fminunc(f,angles,options)

% some plotting at the end
nx1 = radius * cos(angles1) + xCenter; 
ny1 = radius * sin(angles1) + yCenter;
nx2 = radius * cos(angles2) + xCenter; 
ny2 = radius * sin(angles2) + yCenter;
plot(nx1, ny1, 'g*');
plot(nx2, ny2, 'm*');

For other optimization tools type help optim ...

 Optimization Toolbox
 Version 7.2 (R2015a) 09-Feb-2015

 Nonlinear minimization of functions.
   fminbnd      - Scalar bounded nonlinear function minimization.
   ...