How to find the maximum of a 2d function in specific interval?

I have a 2d function generated from loops, now I want to find the maximum value of that function in specific interval ( -1<(x,y)<1). How it is possible in Matlab?

syms x y
P=0;
for i=1:7
    for j=1:7
        P= legendreP(i-1,x)*legendreP(j-1,y)+P;
    end
end

Answer 1

[Dyuman Joshi]

Finding an extremum of a function in a bounded interval requires the use of numerical methods -

syms x y
P = 0;
for i = 1:7
    for j = 1:7
        P = legendreP(i-1, x) * legendreP(j-1, y) + P;
    end
end

In order to maximize P, we will minimize -P

%% Convert the symbolic expression into a anonymous function
%% of a single variable
fun = matlabFunction(-P, "Vars", {[x, y]})

fun = function_handle with value:
    @(in1)in1(:,1).*(-1.1e+1./8.0)-in1(:,2).*(1.1e+1./8.0)-(in1(:,1).*(1.5e+1./8.0)-in1(:,1).^3.*(3.5e+1./4.0)+in1(:,1).^5.*(6.3e+1./8.0)).*(in1(:,2).^2.*(1.05e+2./1.6e+1)-in1(:,2).^4.*(3.15e+2./1.6e+1)+in1(:,2).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,2).*(1.5e+1./8.0)-in1(:,2).^3.*(3.5e+1./4.0)+in1(:,2).^5.*(6.3e+1./8.0)).*(in1(:,1).^2.*(1.05e+2./1.6e+1)-in1(:,1).^4.*(3.15e+2./1.6e+1)+in1(:,1).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,2).^2.*(-1.5e+1./4.0)+in1(:,2).^4.*(3.5e+1./8.0)+3.0./8.0).*(in1(:,1).^2.*(1.05e+2./1.6e+1)-in1(:,1).^4.*(3.15e+2./1.6e+1)+in1(:,1).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,1).^2.*(-1.5e+1./4.0)+in1(:,1).^4.*(3.5e+1./8.0)+3.0./8.0).*(in1(:,2).^2.*(1.05e+2./1.6e+1)-in1(:,2).^4.*(3.15e+2./1.6e+1)+in1(:,2).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,1).*(3.0./2.0)-in1(:,1).^3.*(5.0./2.0)).*(in1(:,2).*(3.0./2.0)-in1(:,2).^3.*(5.0./2.0))+in1(:,2).*(in1(:,1).*(3.0./2.0)-in1(:,1).^3.*(5.0./2.0))+in1(:,1).*(in1(:,2).*(3.0./2.0)-in1(:,2).^3.*(5.0./2.0))-in1(:,1).*in1(:,2)+(in1(:,1).*(3.0./2.0)-in1(:,1).^3.*(5.0./2.0)).*(in1(:,2).^2.*(3.0./2.0)-1.0./2.0)+(in1(:,2).*(3.0./2.0)-in1(:,2).^3.*(5.0./2.0)).*(in1(:,1).^2.*(3.0./2.0)-1.0./2.0)-(in1(:,1).^2.*(1.05e+2./1.6e+1)-in1(:,1).^4.*(3.15e+2./1.6e+1)+in1(:,1).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1).*(in1(:,2).^2.*(1.05e+2./1.6e+1)-in1(:,2).^4.*(3.15e+2./1.6e+1)+in1(:,2).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-in1(:,1).*(in1(:,2).^2.*(3.0./2.0)-1.0./2.0)-in1(:,2).*(in1(:,1).^2.*(3.0./2.0)-1.0./2.0)+(in1(:,2).*(3.0./2.0)-in1(:,2).^3.*(5.0./2.0)).*(in1(:,1).*(1.5e+1./8.0)-in1(:,1).^3.*(3.5e+1./4.0)+in1(:,1).^5.*(6.3e+1./8.0))+(in1(:,1).*(3.0./2.0)-in1(:,1).^3.*(5.0./2.0)).*(in1(:,2).*(1.5e+1./8.0)-in1(:,2).^3.*(3.5e+1./4.0)+in1(:,2).^5.*(6.3e+1./8.0))-in1(:,2).*(in1(:,1).*(1.5e+1./8.0)-in1(:,1).^3.*(3.5e+1./4.0)+in1(:,1).^5.*(6.3e+1./8.0))-in1(:,1).*(in1(:,2).*(1.5e+1./8.0)-in1(:,2).^3.*(3.5e+1./4.0)+in1(:,2).^5.*(6.3e+1./8.0))-(in1(:,1).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,2).^2.*(3.0./2.0)-1.0./2.0)-(in1(:,2).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,1).*(1.5e+1./8.0)-in1(:,1).^3.*(3.5e+1./4.0)+in1(:,1).^5.*(6.3e+1./8.0))+(in1(:,1).*(3.0./2.0)-in1(:,1).^3.*(5.0./2.0)).*(in1(:,2).^2.*(-1.5e+1./4.0)+in1(:,2).^4.*(3.5e+1./8.0)+3.0./8.0)+(in1(:,2).*(3.0./2.0)-in1(:,2).^3.*(5.0./2.0)).*(in1(:,1).^2.*(-1.5e+1./4.0)+in1(:,1).^4.*(3.5e+1./8.0)+3.0./8.0)-(in1(:,1).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,2).*(1.5e+1./8.0)-in1(:,2).^3.*(3.5e+1./4.0)+in1(:,2).^5.*(6.3e+1./8.0))-in1(:,2).*(in1(:,1).^2.*(-1.5e+1./4.0)+in1(:,1).^4.*(3.5e+1./8.0)+3.0./8.0)-in1(:,1).*(in1(:,2).^2.*(-1.5e+1./4.0)+in1(:,2).^4.*(3.5e+1./8.0)+3.0./8.0)-(in1(:,1).*(1.5e+1./8.0)-in1(:,1).^3.*(3.5e+1./4.0)+in1(:,1).^5.*(6.3e+1./8.0)).*(in1(:,2).*(1.5e+1./8.0)-in1(:,2).^3.*(3.5e+1./4.0)+in1(:,2).^5.*(6.3e+1./8.0))-(in1(:,2).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,1).^2.*(-1.5e+1./4.0)+in1(:,1).^4.*(3.5e+1./8.0)+3.0./8.0)-(in1(:,1).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,2).^2.*(-1.5e+1./4.0)+in1(:,2).^4.*(3.5e+1./8.0)+3.0./8.0)-(in1(:,1).*(1.5e+1./8.0)-in1(:,1).^3.*(3.5e+1./4.0)+in1(:,1).^5.*(6.3e+1./8.0)).*(in1(:,2).^2.*(-1.5e+1./4.0)+in1(:,2).^4.*(3.5e+1./8.0)+3.0./8.0)+(in1(:,2).*(3.0./2.0)-in1(:,2).^3.*(5.0./2.0)).*(in1(:,1).^2.*(1.05e+2./1.6e+1)-in1(:,1).^4.*(3.15e+2./1.6e+1)+in1(:,1).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,2).*(1.5e+1./8.0)-in1(:,2).^3.*(3.5e+1./4.0)+in1(:,2).^5.*(6.3e+1./8.0)).*(in1(:,1).^2.*(-1.5e+1./4.0)+in1(:,1).^4.*(3.5e+1./8.0)+3.0./8.0)+(in1(:,1).*(3.0./2.0)-in1(:,1).^3.*(5.0./2.0)).*(in1(:,2).^2.*(1.05e+2./1.6e+1)-in1(:,2).^4.*(3.15e+2./1.6e+1)+in1(:,2).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-in1(:,1).^2.*(6.9e+1./1.6e+1)+in1(:,1).^3.*(2.5e+1./4.0)+in1(:,1).^4.*(2.45e+2./1.6e+1)-in1(:,1).^5.*(6.3e+1./8.0)-in1(:,1).^6.*(2.31e+2./1.6e+1)-in1(:,2).^2.*(6.9e+1./1.6e+1)+in1(:,2).^3.*(2.5e+1./4.0)+in1(:,2).^4.*(2.45e+2./1.6e+1)-in1(:,2).^5.*(6.3e+1./8.0)-in1(:,2).^6.*(2.31e+2./1.6e+1)-in1(:,2).*(in1(:,1).^2.*(1.05e+2./1.6e+1)-in1(:,1).^4.*(3.15e+2./1.6e+1)+in1(:,1).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-in1(:,1).*(in1(:,2).^2.*(1.05e+2./1.6e+1)-in1(:,2).^4.*(3.15e+2./1.6e+1)+in1(:,2).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,2).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,1).^2.*(1.05e+2./1.6e+1)-in1(:,1).^4.*(3.15e+2./1.6e+1)+in1(:,1).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-(in1(:,1).^2.*(-1.5e+1./4.0)+in1(:,1).^4.*(3.5e+1./8.0)+3.0./8.0).*(in1(:,2).^2.*(-1.5e+1./4.0)+in1(:,2).^4.*(3.5e+1./8.0)+3.0./8.0)-(in1(:,1).^2.*(3.0./2.0)-1.0./2.0).*(in1(:,2).^2.*(1.05e+2./1.6e+1)-in1(:,2).^4.*(3.15e+2./1.6e+1)+in1(:,2).^6.*(2.31e+2./1.6e+1)-5.0./1.6e+1)-1.0./8.0

%Bounds on the independent variables
lb = [-1; -1];
ub = [1; 1];

%Initial point to start the numerical process of finding the minima
%via iterative methods
xy = [0 0];

As we are dealing with a constralined non-linear multivariable function, fmincon is the best option. Note that fmincon() requiers Optimization toolbox.

[xy_max, fun_min] = fmincon(fun, xy, [], [], [], [], lb, ub, [])

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

xy_max = 1×2

0.3863 0.3863

fun_min = -1.3263

%The corresponding maximum value
fun_max = -fun_min

fun_max = 1.3263

Answer 2

[recent works]

You can find the maximum value of your 2D function within a specific interval in MATLAB by using the following steps:

1. Define the symbolic variables and the function.
2. Define the interval of interest.
3. Use the fmincon function to find the maximum value within the specified interval.

how you can modify your code to find the maximum value in the interval (-1 < x, y < 1):

syms x y
P = 0;
for i = 1:7
    for j = 1:7
        P = legendreP(i-1, x) * legendreP(j-1, y) + P;
    end
end

% Define the function you want to maximize
f = -P;  % Use negative sign to find the maximum

% Define the constraint for the interval (-1 < x, y < 1)
lb = [-1; -1];  % Lower bound for x and y
ub = [1; 1];    % Upper bound for x and y

% Initial guess (optional)
x0 = [0; 0];  % You can choose a different initial guess if needed

% Solve for the maximum within the interval
options = optimset('Display', 'off');  % Turn off display for fmincon
[x_max, f_max] = fmincon(f, x0, [], [], [], [], lb, ub, [], options);

Error using optimfcnchk
If FUN is a MATLAB object, it must have an feval method.

Error in fmincon (line 435)
funfcn = optimfcnchk(FUN,'fmincon',length(varargin),funValCheck,flags.grad,flags.hess,false,Algorithm);

% Display the result
fprintf('Maximum value: %f\n', -f_max);  % Use negative sign to get the actual maximum
fprintf('Coordinates of the maximum: x = %f, y = %f\n', x_max(1), x_max(2));

• We define the function you want to maximize as f and use the negative sign to find the maximum (since fmincon minimizes by default).
• We define the constraints for the interval (-1 < x, y < 1) by specifying lower and upper bounds lb and ub.
• We use fmincon to find the maximum value within the specified interval.
• Finally, we display the maximum value and the coordinates where it occurs.

[Walter Roberson]

Your f is a symbolic expression. fmincon requires the handle to a function (or, more obscurely, a character vector or string scalar that names a public function -- not a local or nested function.)

Also, fmincon never promises to find a global minima, not even when there are not constraints. fmincon finds a local minima that the initial point is in the "basin of attraction" of.

Answer 3

[Walter Roberson]

syms x y
P=0;
for i=1:7
    for j=1:7
        P= legendreP(i-1,x)*legendreP(j-1,y)+P;
    end
end

Pe = collect(expand(P), [x y])

Pe = 

partialx = solve(diff(Pe,x),x)

partialx = 

Pe2 = collect(expand(subs(Pe, x, partialx)),y)

Pe2 = 

partialy = arrayfun(@(EXPR) double(solve(diff(EXPR))), Pe2, 'uniform', 0)

partialy = 5×1 cell array
    {5×1 double}
    {5×1 double}
    {5×1 double}
    {5×1 double}
    {5×1 double}

You now have a situation where for each of the 5 symbolic partialx values, you have 5 specific double precision partialy values that are critical points for that particular partial x. You would have to substitute the values back into the specific partialx to get the full x value, and then you would evaluate the function at those x and y values; you would put everything together and determine the maximum.

Your P is a multinomial of degree 6 in each of x and y, so there are 5 x 5 = 25 critical points to examine to determine if they are maxima or minima. But you can skip the ones involving imaginary coefficients.

You can seldom find symbolic solutions for the roots of equations of degree 5 or higher, so you are not likely to be able to get to closed form solutions.