Asked by ASR97
on 2 Jul 2019

I am currently working on a project where I need to solve the following equation in order to obtain the lagrangian multiplier, , where I have given vector (where each elemnt is a random number chosen uniformly from [0,1]) and a given value λ:

I am attempting to use numerical methods on MATLAB for the past few days to try and get different values for as I vary the parameter , and plot a graph of ρ against for .

So far I have been obtaining very inconsistent answers and using iterative processes is proving inefficient. This is my latest attempt at coding it this, but I am not having much luck with the function:

phi=[];

LOAD=[];

P=[];

L=[];

mu = rand(1000,1);

MU = sum(mu);

increment = 1/(length(mu));

for j = 0:increment:1

L = j*MU;

LOAD(1/j)= L;

for i = 1:length(mu)

Total(i) = sum((sqrt((mu(i)/(L*phi))))*heaviside(mu(i) - (sqrt((mu(i)/(L*phi))))) + (mu(i)*heaviside((sqrt((mu(i)/(L*phi))) - mu(i)))));

end

syms phi;

P(1/j) = vpasolve(sum(Total)= MU - L, phi)

end

plot(P,LOAD)

This code returns :

The expression to the left of the equals sign is not a valid target for an assignment.

I also previously tried the following, but it waas flat out just giving me incorrect results (I was getting values in the imaginary space)

mu = rand(1000,1);

MU = sum(mu);

increment = 0.01;

for j = 0:increment:1

L = j*MU;

counter = j/increment

if j == 0

X(1) = L/MU;

else

X(round(j/increment + 1)) = L/MU;

end

phi = 0.2;

Total = [];

diff = MU - L;

A = abs(sum(Total) - diff);

B = 100;

while B > 0.2

for i = 1:length(mu)

Total(i) = sum((sqrt((mu(i)/(L*phi))))*heaviside(mu(i) - (sqrt((mu(i)/(L*phi))))) + (mu(i)*heaviside((sqrt((mu(i)/(L*phi))) - mu(i)))));

end

B = abs(sum(Total) - diff);

if B < 0.2

continue

else

if phi - 0.0001 < 0

break

else

phi = phi - 0.0001

end

end

end

for i = 1:length(mu)

summed(i) = sqrt((L*mu(i)*phi) - 1)*heaviside(sqrt(L*mu(i)*phi) - 1);

end

if j == 0

delay(1) = (1/L)*sum(summed);

else

delay(round(j/increment) + 1) = (1/L)*sum(summed);

end

end

plot(X,delay)

If anyone can guide me to a more efficient method of solving this equation or correcting my errors in yntax, that would be greatly appreciated. I have also tried using Mathematica to solve this problem, but it was having issues with the Heaviside Theta Function.

Answer by David Goodmanson
on 3 Jul 2019

Edited by David Goodmanson
on 3 Jul 2019

Hi A7,

The code below finds a solution. Instead of looking for phi0 directly it uses the variable

beta = 1/(lambda*phi0)

Since all the parameters mui,lambda,phi0 are positive,

heaviside( mui -sqrt(mui/(lambda*phi0)) )

= heaviside( mui^2-mui/(lambda*phi0) )

= heaviside( mui-1/(lambda*phi0 )

= heaviside( mui-beta )

Since 0<=mui<=1, the sought-after beta pretty much has to lie in the same range. This narrows down the search range.

You can put a for-loop around this code (functions generally have to go to the end) to obtain phi0 as a function of lambda.

mui = rand(1,10);

mu = sum(mui);

lambda = mu/2 % for example

beta0 = fzero(@(beta) Sfun(beta,mui,mu,lambda), [0 1]);

phi0 = 1/(lambda*beta0);

% check that the solution works with the orignal equation

chek = sum( sqrt(mui/(lambda*phi0)).*heaviside(mui-sqrt(mui/(lambda*phi0))) ...

+ mui.*heaviside(sqrt(mui/(lambda*phi0))-mui)) - (mu - lambda);

max(abs(chek)) % should be small

function S = Sfun(beta,mui,mu,lambda)

S = sum( sqrt(beta*mui).*heaviside(mui-beta) ...

+ mui.*heaviside(beta-mui)) - (mu - lambda);

end

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.