Definite integral with a parameter in the bound

I'm trying to integrate the following function f = ((W-x)^2)*((x^2-1)^(1/2))*x where W is a parameter, from 1 to W with the following code:

syms x W

f = ((W-x)^2)*((x^2-1)^(1/2))*x;

F = int(f,x,1, W)

Now, I don't understand why it doesn't work. Thank you in advance for the help.

[John D'Errico]

Actually, an analytical solution does exist. It just takes a little effort. Sometimes you need to push to get a symbolic solution.

Answer 1

[John D'Errico]

syms X W
F = ((W-X)^2)*((X^2-1)^(1/2))*X;
int(F,X,1, W)

ans = 

As you found, int fails to generate a result. That does not mean a closed form solution does not exist. It does suggest that possibility though. It may be just that int is unable to solve the problem. In fact, I checked, and Wolfram Alpha also fails. Still that does not mean nothing can be done.

expand(F)

ans = 

And there we see the integral can be brokn into three terms. Are they all impossible to solve? Clearly not. The W is a constant, so I'll ignore that for now, just computing indefinite integrals.

int(X*sqrt(X^2 - 1),X)

ans = 

int(X^3*sqrt(X^2 - 1),X)

ans = 

It is an integral of this intermediate general form (with an even power of X outside) that seems to be the sticky part, at least for int.

int(X^2*sqrt(X^2 - 1),X)

ans = 

In fact though, Alpha actually was able to integrate a function of that general form. The complete mess, integrating F was just a bit too complicated for Alpha to digest.

https://www.wolframalpha.com/input?i=int%28x%5E2*sqrt%28x%5E2-1%29%29

A problem is, there is a term in there you won't like. There is an asin(x) term. But if x is GREATER than 1, the arc sine will be complex. In fact, Alpha is actually nice enough to show a plot of the result, and it does indeed have an imaginary part in the indefinite integral. So that may suggest there is no real solution after all.

Is that the case? Well, actually not. The INDEFINTE integral has an imaginary part. But that imaginary part appears to be constant for x above 1. And so the imaginary part will disappear when you form the definite integral, between bounds of 1 and W, as long as W>=1.

Of course, you can be lazy, and just compute a numerical solution.

f = @(x,w) ((w-x).^2).*((x.^2-1).^(1/2)).*x;
intf = @(w) integral(@(x) f(x,w),1,w)

intf = function_handle with value:
    @(w)integral(@(x)f(x,w),1,w)

format long g
intf(2)

ans = 
         0.312068786948633

Anyway, with some effort, we could probably get MATLAB to compute the definite integral in an analytical form, since it appears an analytical solution does exist. Sometimes that requires finding a transformation that solve did not see. The obvious transformation here might be something like x = sin(u), but you will need to be careful, because this will now move into the complex plane, since x is greater than 1. The problem is doable. It just takes a little more effort at this point. Thus, if we make the desired transformation, we will have

x = sin(u)

dx = cos(u)du

int(x^2*sqrt(x^2-1)*dx) = -i*int(sin(u)^2*cos(u)^2*du)

the limits will change of course. Here we will be moving along the imaginary axis for the integral, where the limits of integration become [pi/2,asin(W)]. Remember, when W is greater than 1, asin(W) is complex.

syms u
I2 = -i*int(sin(u)^2*cos(u)^2)

I2 = 

Now we can try to put it all together. I think I've gotten all of my i's dotted at this point.

resultW = int((X^3 + W^2*X)*sqrt(X^2-1),X,[1,W]) - 2*W*(subs(I2,u,asin(W)) - subs(I2,u,pi/2))

resultW = 

Does this result agree with integral, when we did it in a numerical form?

intF = matlabFunction(resultW)

intF = function_handle with value:
    @(W)(W.^2-1.0).^(3.0./2.0).*(W.^2.*4.0+1.0).*(2.0./1.5e+1)-W.*(pi.*6.25e-2i-asin(W).*1.25e-1i+W.*sqrt(-W.^2+1.0).*1.25e-1i-W.^3.*sqrt(-W.^2+1.0).*2.5e-1i).*2.0

intF(1) % should be 0

ans = 
     0

intF(2) % should be 0.3120687...

ans = 
         0.312068786948632

Yes, it does. That took a little effort, but it seeems to be working. As we hoped, all of the imaginary terms drop out.

[Antonino Roccaforte]

I really appreciate your answer thank you! But I've seen that Wolfram apha produces indeed the correct result without much effort in terms of a polynomial, a square root and a logarithm (that is integral_1^W (W - x)^2 sqrt(x^2 - 1) x dx = 1/60 (15 W log(sqrt(W^2 - 1) + W) + sqrt(W^2 - 1) (2 W^4 - 9 W^2 - 8)) for W>=1). I think the problem here in Matlab is to declare from the start that the paramater W must be >1, but i don't know how to declare it since I'm only a beginner with Matlab.

[John D'Errico]

I'm not really a beginner, though not as much of an expert as I want. (I doubt that will ever happen anyway.) But this should be sufficient, at least if int will manage a result.

syms X real
syms W real
assume(W > 1)
F = ((W-X)^2)*((X^2-1)^(1/2))*X;
int(F,X,1, W)

ans = 

And int is still stymied. Honestly, I DID try that before I went to the effort of splitting things apart. :)

The form that Alpha generates is I assume equally valid. The presence of a log is not amazing, as often these things have many ways the same expression can be written.

[Antonino Roccaforte]

Thank you very much for your time! It seems that in this case doing the calculation by hand is the better choice! I've tried with the trasformation x=cosh(t) in the case 1<x<W and t>0 and apparently it works!

Answer 2

[Antonino Roccaforte]

It seems a bit strange, but I've found that a little substitution in the piece that causes the major problem make int to produce the correct result. It is sufficient to put t=x+1;

syms x real

syms W real

F =-2*W*(x-1)^2*sqrt(x^2-2*x);

G=W^2*x*sqrt(x^2-1)+x^3*sqrt(x^2-1);

int(F,x,2,W+1)+int(G,x,1,W)

That's crazy.

[Paul]

I think this is similar, if not the same, letting Matlab do the work for the substitution.

syms X real
syms W real
assume(W > 1)
F = ((W-X)^2)*((X^2-1)^(1/2))*X;
syms U real
Fu = subs(F,X,U-1);
intF(W) = int(Fu,U,2,W+1)

intF(W) = 

intF(1)

ans = 

0

intF(2)

ans = 

vpa(intF(2))

ans = 

0.31206878694863289560703390535281