CODEWITHSUNDEEP
matlabsymbolic-math
Santosh Tiwari
Santosh Tiwari

Reputation: 1216

matlab symbolic integration with real gives a complex answer

I wish to get analytical (closed form) solution to the following integral in Matlab. However, Matlab gives me an answer with real & imaginary parts. How do I get it to produce an answer with just the "real" part. Here is the full code.

close all;
clear all;
clc;

syms t real;
syms thetak real;
syms sik real;
syms tbar real;
syms sjk real;

expr = exp(-thetak*((t-sik)^2 + (t-sjk)^2));
Bijk_raw = int(expr,t,0,1);
Bijk = simplify(collect(expand(Bijk_raw)));
fprintf('Bijk is as follows...\n');
pretty(Bijk);

Upvotes: 0

Views: 180

Answers (2)

Jean Marie Becker
Jean Marie Becker

Reputation: 158

The answer you obtain (if you have a Matlab version similar to mine) and that I reproduce here :

  /               /                     2 \ 
  |  1/2   1/2    |   thetak (sik - sjk)  | 
- | 2    pi    exp| - ------------------- | 
  \               \            2          / 
 /    /  1/2          1/2                 \ 
 |    | 2    (-thetak)    (sik i + sjk i) | 
 | erf| --------------------------------- | i - 
 \    \                 2                 / 

    /  1/2          1/2                       \   \ \ 
    | 2    (-thetak)    (sik i + sjk i - 2 i) |   | | 
 erf| --------------------------------------- | i | | / 
    \                    2                    /   / / 
             1/2 
 (4 (-thetak)   )

gives the impression that you have complex number i everywhere.

But in fact it is a false impression due to (-thetak)^(1/2).

Indeed, taking the square root of a negative number will generate a "i" which in turn will "kill" the other "i"s that are "in contact" with it. This cancellation will occur at different places due to the fact that (-thetak)^(1/2) can be found :

1) inside the erf expressions and

2) as a common denominator (last line).

Verify that rule i^2=-1 applies everywhere leaving no chance to the survival of any "i"...

Finally giving (I have set thetak=s^2 with s>0) :

  /               /                      \ 
  |  1/2   1/2    |   s^2 (sik - sjk)^2   |   
- | 2    pi    exp| - ------------------- | 
  \               \            2          / 
 /    /  1/2                   \ 
 |    | 2    s    (sik  + sjk ) | 
 | erf| ----------------------- |  - 
 \    \           2            / 

    /  1/2                          \   \ \ 
    | 2    s   (sik  + sjk  - 2 )   |   | | 
 erf| ----------------------------- |   | |   /   (4 s)
    \             2                 /   / /

Edit : You could have escaped integration. The idea is to convert the quadratic in $exp(-thetak*((t-sik)^2 + (t-sjk)^2))$ under the so-called "canonical form" which in your case is: $exp(-thetak*(((t-A)^2 + B))/C);$ where $A,B,C$ can be expressed as fonctions of sik and sjk (for example $A=(sik+sjk)/2$); in this way, setting $T=t-A$, you are brought back to a classical Gauss integral with formula :

$$\frac{2}/{\sqrt{\pi}} \int_a^b exp(-t^2} dt ) (erf(b) - erf(a))$$

Upvotes: 1

flawr
flawr

Reputation: 11628

What you get has (up to some constant factors) the form

 1i * (c * erf(1i * a) - erf(c * 1i * (a - 2)))

This consist of two terms of the form

- 1i * erf(1i * x)

Which is also known as the imaginary error function erfi(). It turns out that 

erfi(x) = - 1i * erf(1i * x) = 2/sqrt(pi) * integral(@(t)exp(t.^2),0,x)

So for real values of x your expression is in fact real, and this is the case if thetak >= 0 and sik and sjk are real numbers.

The integral you start out with can be reduced to the integral of exp(-t^2) (using some affine transformation), which famously does not have a "closed form", but it is usually written as 

erf(x) = 2/sqrt(pi) * integral(@(t)exp(-t.^2),0,x)

I highly recommend reading the Wikipedia article on the error function.

Furthermore I recommend using a more beginner friendly CAS than the MATLAB symbolic toolbox. One free and opensource CAS I like to recommend is Maxima.

(This is all written in the MATLAB notation due to the lack of LaTeX here on SO.)

Upvotes: 1

Related Questions