Solve (limit) integral expression (symbolic function) to numerical values in Matlab

Hello everyone,

My expression (after final evaluation) turned out to be a very complicated integral expression. It involves inverse of a function, within function.

I am using symbolic function to solve this. Unfortunately i could not get the answer in terms of expression that I need. Can someone guide me with that ?

Following is the final answer that I am getting. (Note that I need it in terms floating-point format)

    dn=2;
    pfa = .1;
    dref = 1;

    syms x z
    lo_x_lim = erfinv(pfa-1)*sqrt(2);
    hi_x_lim = inf;

    lo_z_lim = 0
    hi_z_lim = dn

    %%Expression

    Q(x,z) = .5 +.5*erf(x/sqrt(2))
    Qinv_ = erfinv(Q - pfa/2) %Note that this need to be changed
    f_step_2 = (Qinv_ - x)/(1-z/dref)
    f_step_3 = Q(f_step_2,0)

    g(x,z) = 2*z/dn^2
    h(x,z) = exp(-x^2/2)

    %combining three equations to form one
    fin = f_step_3*g*h

    % Doing double integral wrt x and z in two steps
    %integrate wrt x
    int_x = int(fin,x,-0.2,hi_x_lim)

    %integrate wrt z

    int_z = int(int_x,z,lo_z_lim,hi_z_lim)

Following is my result: (note that I need it in numeric form!)

int_z =

int(int((z*exp(-x^2/2)*(erf((2^(1/2)*(x - erfinv(erf((2^(1/2)*x)/2)/2 + 9/20)))/(2*(z - 1)))/2 + 1/2))/2, x, 1/5, Inf), z, 0, 2)

Thanks Shan

Answer 1

[Walter Roberson]

In a number of places, you have defined a function of two variables but then reference it as if it were a variable rather than a function. You recognized that in one place, but not the others

    Q(x,z) = .5 +.5*erf(x/sqrt(2))
    Qinv_ = erfinv(Q - pfa/2)

By the second of those two lines, Q is a function but you are not treating it like that, and you are not creating a function as the output. You would need

   Qinv_(x,z) = erfinv(Q(x,z) - pfa/2)

Qinv_, f_step_2, f_step_3, g, h and fin all need to be rewritten in terms of functions with parameters and those parameters need to be passed as appropriate to the lower functions. For example,

fin(x,z) = f_step_3(x, z)*g(x, z)*h(x, z)

I think it unlikely that you will be able to get a closed form expression for either integral; you will probably need to vpa() to get the solution in decimal form.

[Walter Roberson]

You have

f_step_2 = (Qinv_ - x)/(1-z/dref)

where dref = 1, so the denominator is (1-z) . But z ranges from 0 to 2, and when it passes through 1, 1-z would be 0 so you would have a division by 0. Your f_step_2 definition therefore introduces a non-removable singularity, and so the integrals have no defined value and cannot be resolved numerically.