Gaussian Approximation - How to implement this math function in Matlab

Question

I'm having some trouble trying to implement the following math function:

I directly calculated by myself the inverse function of phi(x) for the first equality of the piecewise defined function.

I have the impression that something must have been done wrongly, as the result should be more 'positive' (greater than 0) for each iteration.

I know for sure that this is the exact formula I am supposed to use, so would you please be so kind to give me any feedback about how to solve this?

Many thanks in advance, and BR.

Answer 1

The solution was really simple after all. Glancing through many papers where this algorithm is mentioned, I noticed that the index of both sums don't start in 'j/i =1', but from 'j/i=2', and therefore the exponential is no longer powered by 0.

 l = [0 0.3078 0.27287 0 0 0 0.41933];
 r = [0 0 0 0 0 0.4 0.6];

 sigma = 9.8747;

 mu0 = 2/sigma;

 iterations = 50;

 % Density evolution algorithm depiction for finding the treshold of irregular LDPC codes
 syms x;

 l_idle = zeros(1,length(l));
 r_idle = zeros(1,length(r));

 Q_1 = exp(-0.4527*x^0.86 + 0.0218);
 Q_2 = sqrt(pi/x)*exp((-x/4)*(1-20/(7*x)));

 mv = zeros(1,iterations+1);

 for k=2:length(mv)
    for i = 2:length(l_idle)
         if ((mu0 + (i-1)*mv(k-1)) < 10) 
             l_idle(i) = double(subs(Q_1,x,(mu0 + (i-1)*mv(k-1))));
         else
             l_idle(i) = double(subs(Q_2,x,(mu0 + (i-1)*mv(k-1))));
         end
    end
    lambda = l(2:length(l))*transpose(l_idle(2:length(l_idle)));    

    for j = 2:length(r_idle)
        b = 1-(1-lambda)^(j-1);
        if b < 10, r_idle(j) = subs(0.4527^(-1/0.86)*(0.0218-log(x))^(1/0.86),x,b);
        else,      r_idle(j) = subs(finverse(Q_2,x),x,b);
        end       
    end
    mv(k) = r(2:length(r))*transpose(r_idle(2:length(r_idle)));
end

Many thanks in advance for your support and may you have a nice weekend!