How to find inverse of a self written integral function?

I have defined a function y=iMyF(x) containing an integral written as y=5*log10[4345*\int_0^x du/sqrt(0.3*(1+u)^2+0.7)] so that I can find y for a given value of x. How can I find x for a given value y for this function? I want to use the inverse function to fit a data set using the curve fitting tool.

[Walter Roberson]

I cannot tell from your notation whether

y = 5*log10(4345*(int(sqrt(3/10*((1+u)^2)+7/10), u = 0 .. x)))

or

y = 5*log10(4345*(int(1/sqrt(3/10*((1+u)^2)+7/10), u = 0 .. x)))

?

If it is the second of those then there is the closed form solution,

x = -1+(1/3)*sqrt(3)*sqrt(7)*sinh(arcsinh((1/7)*sqrt(3)*sqrt(7))+(1/43450)*sqrt(3)*10^(1/2+(1/5)*y))

[Rajendra Gupta]

Thanks Walter. Indeed it is the second one. However, the general equation is y = 5*log10(R*(int(1/sqrt(a*((1+u)^2)+(1-a)), u = 0 .. x))). Does it have a closed form solution too? If, it does, that would be great. Thanks again.

[Walter Roberson]

The general solution appears to be

 y = (1/2)*((a^2-a)*exp((1/2)*(-2*exp((1/5)*y*ln(10))*a^(1/2)+R*ln(1/(a*(a^(1/2)+1)^2)))/R)+exp((1/2)*(2*exp((1/5)*y*ln(10))*a^(1/2)+R*ln(a*(a^(1/2)+1)^2))/R)-2*a)/a

[Rajendra Gupta]

Thanks again Walter. I think you mean x= .... (not y= ....).

[Walter Roberson]

Ah, yes, that is correct.

[Rajendra Gupta]

I checked the expression assuming it as x=.... . However, I always get infinity for all values of a & R.

[Walter Roberson]

syms y a R
x(y, a, R) = (exp((a^(1/2)*exp((2592480341699211*y)/5629499534213120) + (R*log(a*(a^(1/2) + 1)^2))/2)/R)/2 - a + (a*exp(-(2*a^(1/2)*exp((2592480341699211*y)/5629499534213120) - R*log(1/(a*(a^(1/2) + 1)^2)))/(2*R))*(a - 1))/2)/a;

fplot( x(y, 3/10, 4345) )

[Rajendra Gupta]

Thanks a million Walter. It works. In fact the other formula you gave also works. It was my mistake that I was careless in defining y. It is y = 5*log10(R*(int(1/sqrt(a*((1+u)^3)+(1-a)), u = 0 .. x))). The difference is in the exponent of (1+u). It is 3 not 2. I wonder if it also has a closed form solution. Thanks again for all the help.

[Walter Roberson]

The integral has a closed form solution in EllipticE for the case x>0, but finding the inverse of it is taking rather some time. Complex numbers are involved.

For the case x < 0, it depends upon where x is relative to a formula involving a that tends towards -infinity as a approaches 0; the parts look like they will be complicated.

[Rajendra Gupta]

In fact x is always greater than zero. So, please just suggest a solution for x>0. Thanks.

[Walter Roberson]

This is turning out to be difficult.

One thing I have found is that for a below approximately 0.28, the integral is imaginary for small values of x. The boundary x gets larger as a gets smaller. It appears to be about

((3*a^(4/3)*(-a+1)^(2/3)+(2*a^(2/3)-3*a^(5/3))*(-a+1)^(1/3)+3*a^2-4*a)*sqrt((-a+1)^(1/3)+a^(1/3))+sqrt(3*(-a+1)^(2/3)*a^(2/3)+(-2*sqrt(3)-3*a+4)*(-a+1)^(1/3)-2*a^(1/3)+3*a^(4/3))*(a+(-a+1)^(1/3)*(1+sqrt(3))*a^(2/3)))/(2*sqrt((-a+1)^(1/3)+a^(1/3))*a)

The equation being solved for is

y = (-15*ln(3)+60*ln(R)-20*ln(a)-10*ln(-a+1)+60*ln(EllipticF(2*3^(1/4)*a^(2/3)*(-a+1)^(1/6)*((-a+1)^(1/3)+a^(1/3))^(1/2)/(a+(-a+1)^(1/3)*(1+3^(1/2))*a^(2/3)), (1/4)*2^(1/2)*(1+3^(1/2)))-EllipticF(2*3^(1/4)*a^(2/3)*(-a+1)^(1/6)*(a^(1/3)*x+(-a+1)^(1/3)+a^(1/3))^(1/2)/((x+1)*a+(-a+1)^(1/3)*(1+3^(1/2))*a^(2/3)), (1/4)*2^(1/2)*(1+3^(1/2)))))/(12*ln(2)+12*ln(5))

[Rajendra Gupta]

That is good information on the limits for a real solution. Am I missing something? I do not see any y in the expression. If it is too difficult to find the analytic inverse then we may have to find it numerically, which we have already done as suggested by Jeff Miller. Thanks again.

[Walter Roberson]

The analytic inverse is difficult, perhaps even impossible (I can't tell yet.) But you could be doing the fzero / fsolve approach on the above equation in EllipticF.

[Benjamin]

This is a very useful topic (inverse of function with integral), however I did not quite understand how it is perfomed in MATLAB. It seems that two approaches can be adopted in MATLAB, one is using finverse command (Walter Roberson answer), and the sencond is what is proposed by Jeff Miller, by using fzero/fsolve.

Anyway, can someone help me to do the inverse of this function, L

% x = 0:0.01:1 (actual data for s)
% x = 0.0125; % (example for s)
G = @(s)  gammainc(5/3, 40*s.^0.6 );
F = @(s)  s - 0.0036*G(s);
L =  integral( @(s) abs(1./sqrt(F(s))) , 0 , x=s , 'ArrayValued',true);
L_inv = ?

your help is very much appreciated.

[Walter Roberson]

Could you confirm that you want the scaled gamma incomplete and not the unscaled?

It looks plausible to me that the sqrt() could be acting on a negative value, leading to complex, but then converted to real by the abs() ? Does that sound accurate?

[Walter Roberson]

I see article https://www.researchgate.net/publication/3388532_Inverse_incomplete_gamma_function_and_its_application . Perhaps some of its ideas could be used as steps in developing the expression you need. It would not be easy though.

Answer 1

[Jeff Miller]

No doubt there are more efficient methods for your particular function, but a quick and dirty general approach is to make a function using fzero. I don't do this often, but I think it would look something like this:

function x = myinverse(y, guess)
  fun = @(z) (iMyF(z)-y);  % parameterized function
  x = fzero(fun,guess)     % guess is your initial pretty good guess for x
end

[Rajendra Gupta]

Thanks Jeff. It works pretty well provided the 'guess' is reasonable. Problem is when y and x are arrays. One then needs to provide 'guess' also as an array, or have a formula included in the function that can provide a pretty good value of guess for each value of y. I would also like to know more efficient methods if at all possible. Nevertheless, for now I am happy as it solved my immediate problem by embedding a formula for guess in the function.