MATLAB: Solving Integral Equation With Unknown Boundary

Question

Hello,

I'm trying to solve the equation below, we know gamma_bar = 100 and K = 0.1973 and we need to find gamma_k which is the lower bound of integral.

I tried to write code by using fsolve and int functions but to be honest, I couldn't managed to run it. In below code, x represents gamma_k and s represents gamma.

Pb = 0.0001;
K = -1.5/log(5*Pb);
gamma_bar = 100;
fx = int((1/K)*((1/x)-(1/s))*(1/gamma_bar)*exp(-s/gamma_bar), s,[x 99999])
x=fsolve(fun,0.01);

As an error message it says "Unrecognized function or variable 'x'" but I'm trying to find x anyway.

Can anyone help me to get gamma_k result ?

Thanks in advance.

Answer 1

Is gammaK both the lower limit of the integral, AND also inside the integral? I would note that you write it differently in the formula you provide versus the code. Of course, the code you provide is not correct anyway, so I'll presume the problem is to solve for gammaK, as the formula indicates.

Next, there is often no good reason to mix symbolic tools and numerical ones. Just because you don't know the value of gammaK does not make it symbolic. This is a common mistake I see made.

But first, I'll try using the symbolic toolbox, since I expect this function is integrable in a useful form.

Note that I'll avoid using the name gamma, as gamma is such a useful tool in MATLAB. Avoid stepping on functions you might need to use later.

% Givens:
Pb = 0.0001;
K = -1.5/log(5*Pb);
gamma_bar = 100;
% unknowns:
syms gammaK real positive
syms gam real positive
kernel = (1/gammaK - 1/gam)/(K*gamma_bar)*exp(-gam/gamma_bar);

Now, don't jump to immediately throw this into a solver. First, is it integrable, or at least is there some solution for the indefinite integral?

int(kernel)
ans =
- (281474976710656*ei(-gam/100))/5554767572841873 - (28147497671065600*exp(-gam/100))/(5554767572841873*gammaK)

Not surprisingly, this is not a difficult problem, so MATLAB finds something, using the special function ei. You might see, I am sneaking up on the solution. testing each part to see if what I have done makes common sense.

As a definite integral, we get something we expect.

gkint = int(kernel,gam,[gammaK,inf])
gkint =
(281474976710656*ei(-gammaK/100))/5554767572841873 + (28147497671065600*exp(-gammaK/100))/(5554767572841873*gammaK)

We can even evaluate this, or turn it into a function in MATLAB.

gkfun = matlabFunction(gkint)
gkfun =
  function_handle with value:
    @(gammaK)ei(gammaK.*(-1.0./1.0e+2)).*5.067268306361389e-2+(exp(gammaK.*(-1.0./1.0e+2)).*5.067268306361389)./gammaK

It appears to cross 1.

gkfun(1)
ans =
       4.8122
       
gkfun(10)
ans =
      0.36613
      
fplot(gkfun,[1,10])
yline(1);

Finally, now that we know how this behaves, we can try to solve for gammaK. We can do that wither symbolically, in the form

vpasolve(gkint == 1,10)
ans =
4.2866521240940762645406270045941

Or we can do it numerically. Use fzero. See that I subtract 1, because fzero looks for a ZERO.

[gksol,fval] = fzero(@(gk) gkfun(gk) - 1,10)
gksol =
          4.28665212409408
fval =
     0

As I said, solve a problem by sneaking up on it. Plot everything. Think about what you get, what you see. Does it make sense? Only at the end should you use a solver to solve for something. Also, you need to understand the difference between symbolic tools and numerical tools. fzero is a numerical tool. It works here only because we used matlabFunction.