David
David

Reputation: 1001

Indefinite integration with Matlab's Symbolic Toolbox - complex solution

I'm using Matlab 2014b. I've tried:

clear all
syms x real
assumeAlso(x>=5)

This returned:

ans =

[ 5 <= x, in(x, 'real')]

Then I tried:

int(sqrt(x^2-25)/x,x)

But this still returned a complex answer:

(x^2 - 25)^(1/2) - log(((x^2 - 25)^(1/2) + 5*i)/x)*5*i

I tried the simplify command, but still a complex answer. Now, this might be fixed in the latest version of Matlab. If so, can people let me know or offer a suggestion for getting the real answer?

The hand-calculated answer is sqrt(x^2-25)-5*asec(x/5)+C.

Upvotes: 1

Views: 271

Answers (2)

horchler
horchler

Reputation: 18484

This behavior is present in R2017b, though when converted to floating point the imaginary components are different.

Why does this occur? This occurs because Matlab's int function returns the full general solution when you ask for the indefinite integral. This solution is valid over the entire domain of of real values, including your restricted domain of x>=5.

With a bit of math you can show that the solution is always real for x>=5 (see complex logarithm). Or you can use more symbolic math via the isAlways function to show this:

syms x real
assume(x>=5)
y = int(sqrt(x^2-25)/x, x)
isAlways(imag(y)==0)

This returns true (logical 1). Unfortunately, Matlab's simplification routines appear to not be able to reduce this expression when assumptions are included. You might also submit this case to The MathWorks as a service request in case they'd consider improving the simplification for this and similar equations.

How can this be "fixed"? If you want to get rid of the zero-valued imaginary part of the solution you can use sym/real:

real(y)

which returns 5*atan2(5, (x^2-25)^(1/2)) + (x^2-25)^(1/2).

Also, as @SardarUsama points out, when the full solution is converted to floating point (or variable precision) there will sometimes numeric imprecision when converting from exact symbolic form. Using the symbolic real form above should avoid this.

Upvotes: 2

Sardar Usama
Sardar Usama

Reputation: 19689

The answer is not really complex.

Take a look at this:

clear all; %To clear the conditions of x as real and >=5 (simple clear doesn't clear that)
syms x;
y = int(sqrt(x^2-25)/x, x)

which, as we know, gives:

y =

(x^2 - 25)^(1/2) - log(((x^2 - 25)^(1/2) + 5i)/x)*5i

Now put some real values of x≥5 to check what result it gives:

n = 1004;             %We'll be putting 1000 values of x in y from 5 to 1004
yk = zeros(1000,1);   %Preallocation
for k=5:n           
    yk(k-4) = subs(y,x,k);    %Putting the value of x
end

Now let's check the imaginary part of the result we have:

>> imag(yk)

ans =

   1.0e-70 *

                   0
                   0
                   0
                   0
   0.028298997121333
   0.028298997121333
   0.028298997121333

%and so on...

Notice the multiplier 1e-70.

Let's check the maximum value of imaginary part in yk.

>> max(imag(yk))

ans =

     1.131959884853339e-71    

This implies that the imaginary part is extremely small and it is not a considerable amount to be worried about. Ideally it may be zero and it's coming due to imprecise calculations. Hence, it is safe to call your result real.

Upvotes: 2

Related Questions