matlab not displaying an approximated value of an improper integral?

I want to calculate an approximate value of an improper integral:

>> fun= 4./(x.*(1-x).*(2-x).*(pi.^2+(log(x./(1-x)).^2))

>> vpa(int(fun,0,1))

then matlab displays the result:

ans =

vpaintegral(4/(x*(log(-x/(x - 1))^2 + 2778046668940015/281474976710656)*(x - 1)*(x - 2)), x, 0, 1).

how can I get the approximate value of this improper integral? thanks in advance.

[Nathan Hardenberg]

Hey Lahcen, I can not reproduce your results. Your fun variable is not correct. It seems there is a bracket missing at the end. But if I add it I do not get the answer you got. Have you done something beforehand? Can you show how you defined x and p? And make sure to "clear 'all'" at the start of the script

clear 'all'
syms x p
fun= 4./(x.*(1-x).*(2-x).*(p.^2+(log(x./(1-x)).^2)))  % added bracket at the end

fun = 

vpa(int(fun,0,1))

ans = 

[Dyuman Joshi]

"But if I add it I do not get the answer you got"

If you check the term given by OP as the result and the final output of your code, they are the same (with a value substituted for p)

@Lahcen

You can try solving numerically, but as there are singularities at x==0 and x==1, the result will not be accurate -

%p^2 value taken from above
p2=2778046668940015/281474976710656;
fun= @(x) 4./(x.*(1-x).*(2-x).*(p2+(log(x./(1-x)).^2)));
integral(fun,0,1)

Warning: Minimum step size reached near x = 1. There may be a singularity, or the tolerances may be too tight for this problem.

ans = 2.7246

[Nathan Hardenberg]

@Dyuman Joshi good catch. I saw that but did not realise that p is probably just this numerical value then

p = 2778046668940015/281474976710656

p = 9.8696

[Lahcen]

I'm sorry, it's pi not p.

[Dyuman Joshi]

I had an inkling because of the value but wasn't sure due to the format of the value -

sqrt(2778046668940015/281474976710656)

ans = 3.1416

[Lahcen]

Thank you al, Mr Nathan Hardenberg, Mr Dyuman Joshi and Mr Steven Lord,  for your helpful responses in addressing my improper integral calculation issue on the Matlab forum. Your assistance was greatly appreciated!. Indeed, I worked by the idea of Mr Nathan Hardenberg.

I'd also like to note that the integral in question is convergent, whereas the numerical methods provided by Matlab seem to be getting stuck without yielding the desired result. It appears that there is a need for the development of more reliable numerical techniques.

Many thanks.

Answer 1

[Nathan Hardenberg]

Here are some ideas, but please still regard my comment.

1. You have two variables x and p. You can not approximate numerically then.
2. "To approximate integrals directly, use vpaintegral instead of vpa. The vpaintegral function is faster and provides control over integration tolerances." [source]
3. Specify that you want to use x as your variable (int(fun, x, [a b]))
4. You can not numerically approximate if your denomenator gets 0 at some point (see example below)

-- EDIT --

Your problem is nr. 4. Your denomenator gets 0 to fix this you can simply choose values close to 0 and 1. See example below:

syms x
p = pi; % edited to be p = pi
fun = 4./(x.*(1-x).*(2-x).*(p.^2+(log(x./(1-x)).^2)))

fun = 

vpa(int(fun, 0.001, 0.999))  % choosing values slightly above 0 and below 1

ans = 

2.0701395395280001703925025323768

But maybe this is not good enough, since the solution can get quite a bit different when getting closer to the limits:

vpa(int(fun, 1e-110, 1 - 1e-13))  % choosing tighter values

ans = 

2.7443512001078913067737891865659

-- EDIT END --

% Example of nr. 4 from above
syms x
f = 1/x;
Fvpaint = vpaintegral(f,x,[0 1])

Error using sym/vpaintegral
Failed precision goal. Try using 'MaxFunctionCalls'.

[Steven Lord]

Let's look at how your function behaves over the limits of integration. I had to add one extra closing parenthesis, which I did at the end. I'm not sure if that's the right place for it.

syms x
fun= 4./(x.*(1-x).*(2-x).*(pi.^2+(log(x./(1-x)).^2)))

fun = 

fplot(fun, [0 1])

That certainly looks like a singularity at x = 1.

try
    vpa(subs(fun, x, 1))
catch ME
    fprintf('subs threw an error:\n%s', ME.message)
end

subs threw an error:
Division by zero.

The value at x = 0 also looks like it's potentially very, very large (if not infinite.)

try
    vpa(subs(fun, x, 0))
catch ME
    fprintf('subs threw an error:\n%s', ME.message)
end

subs threw an error:
Division by zero.

What would I get if I tried to numerically integrate the function?

f = matlabFunction(fun)

f = function_handle with value:
    @(x)4.0./(x.*(log(-x./(x-1.0)).^2+9.869604401089358).*(x-1.0).*(x-2.0))

integral(f, 0, 1)

Warning: Minimum step size reached near x = 1. There may be a singularity, or the tolerances may be too tight for this problem.

ans = 2.7246

What do you think the value of this integral should be and why?

[Lahcen]

It is convergent, according to Bertrand's criterion.

[Dyuman Joshi]

@Lahcen Can you show the calculations you used to determine that the integral in question is convergent via Bertrand's criterion?

[Lahcen]

@Dyuman Joshi: Dear Dyuman, I apologize for the delay in response.

Recall that:

By Bertrand, let a<0, converge if and only if or ()

Now, in the vicinity of 0, we have

, which converges according to Bertrand's criterion.

For 1, it suffices to set x=1-t to find

, which converges according to Bertrand's criterion.

whence the integral is convergent.

[Lahcen]

@Nathan Hardenberg

Thank you very much for the idea of the approximation, I based on your idea. Indeed, I noticed that the approximation via matlab in the vicinity of 0 does not present a problem while in the vicinity of 1 Matlab only approached 1 approximately , so I split the integral into two parts as follows:

then I made a change of variable for the second integral t=1-x, we get

Now we will apply your idea of approximation in the vicinity of 0 for the two integrals via matlab, we find:

[Nathan Hardenberg]

@Lahcen Glad I I could help. Nice idea to use the better numerical precision near 0 to get a better approximation

Now I also got interested in the convergance 😀 I did not find the formula by Bertrand you provided. I only found the following. Always converging when

• (source)
• (source)

But this is obiously somthing else, regarding the area after to infinity. Would be greatly appreciated if you would provide a resource, if you can. But no worries if you can't

[Lahcen]

Hello Dear @Nathan Hardenberg

Attached you will find a reference about A, concerning your question, (but they are references in French).

Source 1

Source 2

Source 3

But using what you found like, one can find what my published. Indeed, Let a<1 and let put

by setting we find,

then

Now, this integral is converge if and only if

which is equivalent .

I hope you like the proof.

[Nathan Hardenberg]

@Lahcen Thanks a lot :D