How can I get the correct answer for int(1/x) with symbolic integration?

I am trying to solve some ODEs using DSOLVE, and while doing this I realised that matlab is consistently losing half the solutions because it is integrating 1/x wrongly. I'll illustrate this here with the simplest possible example:

syms x
fun = 1/x;
int(fun)

ans = 

which is correct, but only for x>0. The correct answer for all values of x (with the exception of 0) should be log(abs(x)). How can I force matlab to use the full solution?

I also tried

syms x real
fun = 1/x;
int(fun)

ans = 

and

syms x
fun = 1/x;
int(fun,IgnoreAnalyticConstraints=true)

ans = 

[Sam Chak]

It is true that

syms x real

%% function
f   = log(abs(x))

f = 

%% derivative of the function w.r.t. x
df  = diff(f)

df = 

%% simplification done
df  = simplify(df, 'steps', 100)

df = 

However, when the absolute value operator is performed on the input x, there is information loss. Specifically, the sign of the value is lost, meaning you can no longer determine whether the original value was positive or negative.

[Torsten]

@Sam Chak

OP thinks that log(abs(x)) is a more general answer as antiderivative to 1/x because it's valid for all x~=0, not only x > 0. But as said: log(x) is also defined for x < 0 (in the complex plane) with derivative 1/x. So log(abs(x)) as antiderivative to 1/x is not more general than log(x) if we allow complex results.

[Georg]

Many thanks for all the attempts to explain this.

However, I will illustrate the problem with a slightly more elaborate example:

clear all
syms y(x)
eqn = diff(y,x) == (2*x+y(x)+2)/(2*x+y(x)-4); 
solutions = dsolve(eqn,Implicit=true)

solutions = 

which is correct, apart from the fact that the logarithms in the first solution miss the absolute value, which can be seen when plotting the solution

C1val = -20:20;

% The solution from matlab
figure
xlim([-4,4]);
ylim([-2,8]);
pbaspect([1 1.25 1])
grid on
syms u
hold on
for i = 1:length(C1val) 
    impl = subs(lhs(solutions(1))-rhs(solutions(1)),[sym("C1"),y],[C1val(i),u]);
    impl = matlabFunction(impl,'Vars',[x,u]);
    fimplicit(impl,"Color","black")
end
hold off

and it should correctly be

figure
xlim([-4,4]);
ylim([-2,8]);
pbaspect([1 1.25 1])
grid on
hold on
% The solution using int(1/x)=log(abs(x))
for i = 1:length(C1val) 
    f2=@(x,u)u-x-log(abs(+u+x.*2.0-2.0)).*2.0+2.0 - C1val(i);
    fimplicit(f2,"Color","black")
end
hold off
xlim([-4,4]);
ylim([-2,8]);
pbaspect([1 1.25 1])
grid on

ANd the root cause seems to me to be the incmplete solution for int(1/x).

Cheers, Georg

[Torsten]

Try

C1val = -20:20;
C1val = C1val - 2*pi*1i;

[Georg]

Hmm, yes, that works, but what is the argument for adding a complex constant to the integration constant of a real indefinite integral? And why exactly -2pi i? Seems a much less intuitive answer than int(1/x)dx = ln |x|. But as I said, it works. Thank you.

Cheers, Georg

[Walter Roberson]

When you take int(1/x, x, -b, -a) assuming a>0, b>0, then if you get out log(-a)-log(-b) then that is (log(a)+2pi*i) - (log(b)+2*pi*i) which is log(a) - log(b) because the 2*pi*i cancel out.

When you take int(1/x, x, -b, -a) assuming a>0, b>0, then if you get out log(abs(-a))-log(abs(-b)) then that is log(a)) - log(b)

When you take int(1/x, b, a) assuming a>0, b>0 then if you might get out log(a)-log(b) . When you take int(1/x, b, a) then you might get out log(abs(a)) - log(abs(b)) but because a > 0 and b > 0 then that is log(a)-log(b) which is the same.

So whether you get log(x) or log(abs(x)) over b to a, then you get the same result provided that a and b are both less than 0 or a and b are both greater than 0.

You only get a different result between log(x) and log(abs(x)) in the case where one of the values is less than zero and the other value is greater than 0. In which case you end up with either a +2pi*i or -2pi*i for the log(x) form and with no +/-2pi*i for the log(abs(x)) form.

Now consider int(1/x, x, -a, a) then if the integral is log(x) you get log(a) - log(-a) which is log(a) - (log(a) + 2pi*i) which is 2pi*i . If the integral is log(abs(x)) then you get log(abs(a)) - log(abs(-a)) which is log(a)-log(a) which is 0.

syms x real
syms a b real
int(1/x, -a, b)

ans = 

limit(ans, b, a)

ans = 

integral(@(x)1./x, -1, 1)

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 3.7e+00. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

ans = 2.3093e-14

fplot(1/x, [-1 1])

for every -x, x>0, there is a corresponding +x such that 1/(-x) = -1/x cancels out the 1/x . However that still leaves x = 0, for which there is no corresponding match, so the exact integral of 1/x over -a to +a has to be inf or nan -- which is not 0 as would be predicted from log(abs(x))

[Torsten]

@Georg

Hmm, yes, that works, but what is the argument for adding a complex constant to the integration constant of a real indefinite integral? And why exactly -2pi i?

I already answered this:

log(x) = log(-x) + 1i*pi  (x < 0)

Here:

-2*log(2*x+y(x)-2) = -2*log(-(2*x+y(x)-2)) - 2*pi*1i

Thus if you subtract 2*pi*1i from C1, you consider the case that 2*x+y(x)-2 is negative in your implicit plot.

To illustrate this, the following example shows how MATLAB treats cases "left to the line y = -2*x + 2" (e.g. if you impose the condition y(0) = 0). As you can see from the output, MATLAB replaces 2*log(-(2*x+y(x)-2)) by 2*log(2*x+y(x)-2) - 2*pi*i. This circumvents the necessity to work with 2*log(abs(2*x+y(x)-2)).

clear all
syms y(x)
eqn = diff(y,x) == (2*x+y(x)+2)/(2*x+y(x)-4); 
cond = y(0)==0;
solution = dsolve(eqn,cond,Implicit=true)

solution = 

syms u
impl = subs(lhs(solution)-rhs(solution),y,u);
impl = matlabFunction(impl,'Vars',[x,u])

impl = function_handle with value:
    @(x,u)u-x+pi.*2.0i-log(u+x.*2.0-2.0).*2.0+log(2.0).*2.0

hold on
fimplicit(impl)
fimplicit(@(x,u)2*x+u-2)
hold off
grid on

Answer 1

[David Goodmanson]

Hello Georg, I agree that log(|x|) is the correct way to do this, consistent with what tables of integrals have to say. (With implicit equations the situation is more complicated but the principle is the same). When 0<a<b,

Int{a,b} (1/x) dx = log(|b|) - log(|a|)

and as long as the Principal Value of the integral is used, this result holds for the other five orderings of a,b,0, including when a and b are on opposite sides of the origin. For a<0<b the PV integral is

      lim eps-> 0  ( Int{a,-eps} (1/x) dx + Int{eps, b} (1/x) dx )

substituting y = -x in the first integral leads to

      log(eps)-log(-a)  +  log(b)-log(eps) = log(|b|) - log(|a|)

(In this case the limit eps->0 is not really needed, but when integrating something like g(x)/x where g is well behaved across the origin, g can be expanded in a power series, g0 + g1*x + ... . Then eps-> 0 makes a lot of terms drop out.)

The results above are an all-real calculation. For the complex plane aspect, consider a horizontal line displaced above the x axis by a small amout i *delta. With

      z = |z|e^(i*theta)         log(z) = log|z| + i*theta,

then as you move along that line from positive x to negative x, you pick up a term of i*pi as you sweep past the origin. So

      log(z) =  log(|z|) + i*pi                        z<0               (1)

and there is an imaginary term hanging around. However, let b<0<a and integrate to the left from a to b. The path has to be continuous in the complex plane so the path of integration can be path 1,

     Int{a,eps} (1/z) dz + Int{C} (1/z) dz + Int{-eps, b} (1/z) dz
                                       _     
                          b----<---eps/.\eps----<----a     

where C is a semicircle of radius eps. But the integral of (1/z), counterclockwise around that semicircle, is i*pi (consistent with (1)). The other two terms are the PV integral. The result is

log(|b|) - log(|a|) + i*pi       Int(path 1).

So what the principal value is saying is: ignore the semicircle, and obtain the all-real solution.

A convenient version for contour integration is obained by considering path 2 that, rather than going over the top with a semicircle, goes underneath with a semicircle. The result is

log(|b|) - log(|a|) - i*pi      Int(path 2)

and

principal value = (1/2) (Int(path 1) + Int(path 2))     

Answer 2

[Torsten]

If you accept complex constants of integration, log(x) as antiderivative to 1/x is correct also for x < 0 since log(x) = log(-x) + 1i*pi.

Answer 3

[Paul]

Seems like the OP might be missing a sign(x), though I'm not really sure how to intrepret this result.

syms x real
d = diff(log(abs(x)))

d = 

int(d)

ans = 

[Torsten]

x > 0: sign(x)/abs(x) = +1/x = 1/x

x < 0: sign(x)/abs(x) = -1/(-x) = 1/x

syms x real
d = simplify(diff(log(abs(x))),'Steps',40)

d = 

[David Goodmanson]

Georg has a good point. For

w = Int{a to b} (1/x) dx,

if neither a or b are zero and the principal value is taken when integrating across the origin, then for a<b<0 and the five other possible ways of ordering a,b,0

w = log(|b|) -log(|a|) 

works for all of them. Books on integration list the indefinite integral as log(|x|) and leave any complex constants of integration such as i*pi as just that, constants of integration.