Evaluating an integral with symbolic integral limits

Hi,

This is my code:

syms s y0 

result = int( ((1+s)/(((1+s)^3-(1+y0)^3)/3 - ((1+s)^2-(1+y0)^2)/2)^0.5), s, y0, 1);

r = solve(result == 20^0.5);

I was trying to evaluate the following integral which has a lower integral limit as a symbolic variable y0. I wanted to evaluate the result of the integration (which I am assuming should be an expression in terms of y0 according to my code) by using the 'solve' function in matlab by equating the equation to be equal to 20^0.5 (Hence, giving a value for y0).

This is what i get for the result:

 result: int((s + 1)/((s + 1)^3/3 - (s + 1)^2/2 + (y0 + 1)^2/2 - (y0 + 1)^3/3)^(1/2), s == y0..1)

But for the last line of the code, I get this error:

Error using mupadengine/feval (line 157)

MuPAD error: Error: The second argument must be of form x or x = a..b. [int]

Error in solve (line 170)
  sol = eng.feval('symobj::solvefull',eqns,vars);

Error in HW3_graph_code (line 7)

r = solve(result == 20^0.5);

Should I use fsolve to solve the relation?

Thank you

Ushnik

[modem sudhakar]

Hi mukherjee,

how you solved this problem. I am facing similar problem as below. Please sugget the way to solve it.

syms x V; Pr=solve( int(x.^2.* exp(-x-1./(x+V)),x,1,inf )-0.5 ==0,V); Warning: Explicit integral could not be found. Error using mupadengine/feval (line 157) MuPAD error: Error: The second argument must be of form x or x = a..b. [int]
Error in solve (line 160) sol = eng.feval('symobj::solvefull',eqns,vars);

Answer 1

[Roger Stafford]

It is not surprising that 'int' was unable to obtain an explicit expression for this integral, which I believe involves incomplete elliptic integrals of the first and/or second kinds. To obtain a numerical solution to your equation, it would be necessary to create a function which can numerically evaluate the difference between your integral and 20^0.5 as a function of y0. Then you can hand the problem to matlab's 'fzero'. Of course you will need to provide an initial estimate for y0, but you can do that by first using the 'plot' function to show roughly where your function crosses zero. Because a complete numerical integration is required for each step in fzero's iterative procedure, your fzero code (and your plot code, too) will tend to be slow.

[Ushnik Mukherjee]

Hi Roger,

So what you mean is that I should try to get a function that does not have any integrals in it and plot it to find where it crosses zero?

Thank you

Ushnik

[Ushnik Mukherjee]

Hi Roger,

I was just wondering if I could numerically evaluate the integral in order to find y0? This is mu function:

function fun1 = exHW3(x,y0)

fun1 = ((1+x)/(((1+x).^3-(1+y0).^3)/3 - ((1+x).^2-(1+y0).^2)/2).^0.5);

end

I want to calculate for what value of y0 will the above expression when integrated with respect to x (with lower limit of y0 and upper limit of 1) will be equal to 20^0.5. Is there an iterative method you can suggest?

Thank you for the help!

Ushnik

[Roger Stafford]

"So what you mean is that I should try to get a function that does not have any integrals in it". No, what I meant was to prepare a function which you could call, 'myintegral', that accepts y0 as a scalar input and proceeds to use it to perform a numerical integration of your function

   ((1+s)/(((1+s)^3-(1+y0)^3)/3 - ((1+s)^2-(1+y0)^2)/2)^0.5)

with respect to s from s = y0 to s = 1 and then subtracts sqrt(20) from this integral afterwards. You can call on one of the quadrature routines such as 'integral' from inside 'myintegral' to do this.

"Is there an iterative method you can suggest?" Yes, I would suggest matlab's 'fzero'. You could have 'fzero' call on 'myintegral' to repeatedly perform integrations while adjusting y0 until a value of 0 is obtained from 'myintegral' - that is, until the difference between the integral and sqrt(20) becomes zero.

Since the variable y0 is not only an integration limit but is also used within the integrand, it will be necessary to repeat the full integration from beginning to end for each new value of y0. That is what will make this procedure rather slow.

[Walter Roberson]

No, fun1 has imaginary components for nearly all points, and if it has any real values at all they would occur by exactly balancing multiple imaginary components. For example one of the sub-expressions is inherently imaginary except for at y0 = -1, but another sub-expression is inherently imaginary outside y0 from 0 to 1.

There appears to be a closed form to the integral involving complex sign, signum, sin, cos, square roots, absolute value, and short polynomials in y0 up to exponent 3. But finding the points at which that closed form is real-valued is a mess.

[Roger Stafford]

No, I definitely disagree with you, Walter! As long as y0 is kept greater than zero, the quantity inside the square root in the integrand will always be non-negative and will not therefore produce imaginary values. This quantity can be simplified:

((1+s)^3-(1+y0)^3)/3 - ((1+s)^2-(1+y0)^2)/2 =

(s^2-y0^2)/2+(s^3-y0^3)/3

which is obviously non-negative provided s >= y0, which will hold true over the entire stated integration interval.

The integral clearly diverges to plus infinity as y0 approaches 0 and converges to zero as y0 approaches 1, so at least at one point for y0 between 0 and 1 it must cross the value sqrt(20) (it must in fact cross every possible positive value,) and the problem must therefore have a real solution.

[Walter Roberson]

This is an odd situation, in that the general form of the integral, involving EllipticE and EllipticF functions, has EllipticF calls that have imaginary arguments starting at y0 = 1/2 and the situation gets even more complicated starting at y0 = 2/3. And yet, now that I look at what you are pointing out, I agree that only positive quantities are being integrated under the condition you note.

I find that the round-off error with N digits of calculation leaves an imaginary residue on the order of 10^(-N/2) so with 16 digits an imaginary quantity on the order of 10^(-8) is left. That was throwing off my efforts.

[Roger Stafford]

@Walter: I am guessing that the elliptic integral expressions you obtained have imaginary components that end up canceling out in the definite integrals for positive y0 to leave only real values. That is one of the difficulties that turn up in complicated symbolic solutions. A simple situation is sometimes obscured by the necessity of producing very general expressions without the benefit of human intervention.

@Usnhik: As to the numerical difficulties Walter noted, that suggests to me that you really ought to make a change of variable in performing your integration. In place of s, the variable t should be defined as:

 t = sqrt(s-y0)

Doing this will have the effect of eliminating the square root singularity that occurs at the lower limit of integration where s starts out equal to y0. This should lend itself to more accurate numerical integration results.

Also it would be helpful if the integrand expression were simplified along the lines I have previously mentioned:

 ((1+s)^3-(1+y0)^3)/3 - ((1+s)^2-(1+y0)^2)/2 =
 (s^2-y0^2)/2+(s^3-y0^3)/3 =
 (s-y0)*((s+y0)/2+(s^2+s*y0+y0^2)/3)

The (s-y0) factor will disappear when the above variable t is introduced and that will remove the singularity at s = y0, allowing you to obtain better numerical accuracy.

[Ushnik Mukherjee]

Hi Roger,

For some reason I did not get any notification of your posts! Itried what you suggested and simplified the expression and did a change of variables.

clc;

clear all;

clear variables;

close all;

y0 = 0:0.001:1;

a = 20^0.5;

for i= 1:length(y0)

      fun2 = @(t) ((2*(1 + t.^2 + y0))./((t.^2 + 2.*y0)/2 + (2.*y0.^2 + t.^4 + 
  2.*y0.*(t.^2))/3).^0.5);

      q = integral(fun2, 0, (1-y0(i))^0.5);

      if q-a< 0.00001
          y_initial = y0(i)
          i

          break

      end

end

I keep getting this error when i try to run the code!

Error using  * 
Inner matrix dimensions must agree.

Error in @(t)((2*(1+t.^2+y0))./((t.^2+2.*y0)/2+(2.*y0.^2+t.^4+2.*y0*(t.^2))/3).^0.5)

Error in integralCalc/iterateScalarValued (line 315)
                fx = FUN(t);

Error in integralCalc/vadapt (line 133)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 76)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in HW3_y0_calc (line 30)
    q = integral(fun2, 0, (1-y0(i))^0.5);

Where am I going wrong?

I really appreciate your help in this! Ushnik

[Ushnik Mukherjee]

I used it ('ArrayValued', true) then it does not give me an error, but the loop just keeps on going without terminating!

[Ushnik Mukherjee]

Hi Roger,

I got the error sorted out! Thank you very much for all your help! Just one last question, I tried using your substitution of a new variable, and I also evaluated the integral without the substitution. The solutions are a bit different, why is that?

Regards

Ushnik

[modem sudhakar]

Hi mukherjee,

how you solve this problem. I am facing similar problem as below. Please sugget the way to solve this problem.

syms x V; Pr=solve( int(x.^2.* exp(-x-1./(x+V)),x,1,inf )-0.5 ==0,V); Warning: Explicit integral could not be found. Error using mupadengine/feval (line 157) MuPAD error: Error: The second argument must be of form x or x = a..b. [int]

Error in solve (line 160) sol = eng.feval('symobj::solvefull',eqns,vars);