Максимизация функции в Matlab (использование fsolve и производных)

Вот пример того, как это сделать с помощью символического набора инструментов:

%// (there's probably a way to generalize this as well)
dAdp = cellstr(char(...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p1')) ';'],...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p2')) ';'],...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p3')) ';']))

%// convert to proper vector equation
dAdp = regexprep(dAdp, 'p([0-9])', 'p\($1\)');

%// convert to function handle
F = str2func( strcat('@(p) [', dAdp{:}, ']') );

но я бы не рекомендовал делать это таким образом (это совершенно нечитаемо и подвержено ошибкам).

Вы можете написать правильную функцию и оценить dAdp выше, используя символический набор инструментов, но я также не рекомендую делать это таким образом (это просто очень медленно).

Я сторонник вычислений, потому что доверяю компьютерам только то, для чего они созданы: расчеты. Выводы я делаю на бумаге, за исключением чрезмерно простых, но длинных и нудных (можно возразить, что это такой случай, мне до сих пор нравится делать это на бумаге, потому что я также люблю упражнять свои мозг :).

Я предлагаю вам сделать это и/или постоянно проверять себя с помощью символического набора инструментов. ИМХО, вы должны использовать его как помощь, а не как основной двигатель.

Итак, поехали. Вот снова ваша функция, на этот раз в другой форме, с немного большим количеством мозгов:

function A = myfun(p)    
    R = 0.1; 
    u = [500; 400; 300];     
    C = sum(1./u);
    B = C * sum(p) * (1 + R*C*sum(p.^2));

Итак, вы хотите решить ∇A(p) = 0. Лучший способ - применить больше мозгов. Вы можете убедиться, что векторная производная равна:

function F = Aprime(p)
    R = 0.1;
    u = [500; 400; 300];
    C = sum(1./u);
    F = C*( C*R*sum(p.^2) + 2*C*R*p*sum(p) + 1 );

который вы можете решить с большим количеством мозгов:

C·( CR·Σp² + 2CR·p·Σp + 1 ) = 0
                Σp² + 2p·Σp = -1/(CR)

вектор == скаляр: это означает, что все элементы в p равны.
Подставляя q = p1 = p2 = p3, тогда

 3q² + 2q·3q = -1/(CR)
         9q² = -1/(CR)

⇒ q = ⅓·√(-1/(CR))

(что указывает на проблему), или в fsolve вот так:

fsolve(@Aprime, [1 1 1])

что немедленно приведет к

fsolve завершен, потому что вектор значений функции в начальной точке близок к нулю, если измерять значение допуска функции по умолчанию, и проблема кажется регулярной, если измерять градиент.

что действительно указывает на беду.

Поскольку вы указали, что вас интересует максимум, а функция, по-видимому, не имеет стационарных точек, единственный вывод состоит в том, что функция не имеет ни максимума, ни минимума, если только вы не наложите ограничения на p. Действительно, если уменьшить размерность до 2 и построить график:

R = 0.1;
u = [500; 400; 300];
C = sum(1./u);
B = @(p) C * sum(p) .* (1 + R*C*sum(p.^2));

[p1,p2] = meshgrid(-10:0.1:10);
surf(p1,p2,reshape(B([p1(:) p2(:)].'), size(p1)), 'edgecolor', 'none')

... экстремума нет.