Errors when solving equation numerically with a Heaviside Theta

Question

ASR97 2019-7-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/469933-errors-when-solving-equation-numerically-with-a-heaviside-theta

编辑： David Goodmanson 2019-7-3

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2019-7-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/469933-errors-when-solving-equation-numerically-with-a-heaviside-theta#answer_381726

编辑：David Goodmanson 2019-7-3

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