## Errors when solving equation numerically with a Heaviside Theta

### ASR97 (view profile)

on 2 Jul 2019
Latest activity Edited by David Goodmanson

### David Goodmanson (view profile)

on 3 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=[];
P=[];
L=[];
mu = rand(1000,1);
MU = sum(mu);
increment = 1/(length(mu));
for j = 0:increment:1
L = j*MU;
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
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.

R2017b

### David Goodmanson (view profile)

on 3 Jul 2019
Edited by David Goodmanson

### David Goodmanson (view profile)

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