Integration, part 1

Mathcalc8 has three functions for integration. They are called 'int', 'simp' and 'integrate'. Since the function between 'int' and 'integrate' is the same except that 'int' can evaluate complex function whereas 'integrate' can only evaluate real function, we prefer to use 'int' and 'simp' only. In the following post, we will show how these two functions are used by calculating a lot of integrals.

To find the answer for an integral of a function, it is the same as calculating the area that is below the curve stretched by the function. For example, a circle's formula is x² + y² = 1. Suppose we only consider the first quadrant (i.e. region of positive x and y), then the area under the curve of circle is 1/4 of the area of the whole circle :

Input 1 : Pi/4
Ans. 0.78539816
 

To calculate this area, we can use the following integral formula :

∫01 √(1 - x2) dx

The simplest method to calculate this integral is to cut the region inside the circle into many narrow rectangular strips and then sum up all the area of those strips. Below is a formula to implement this algorithm by dividing the region to 50 strips :

Input 2 : f='sqrt(1-x*x)'
 

Input 3 : n=50
 

Input 4 : sum(table('f',id(n)/n))*(1/n)
Ans. 0.79456713
 

The factor 1/n is the width of the strip. Although the answer is only approximately equal to the exact value, this simple algorithm provides a direct way to evaluate any definite integral. (In many textbook, this is actually the definition for definite integration when n → ∞ .)

Next, we try to use the function 'simp' on the same problem :

Input 5 : simp('f',0,1,0.02)
Ans. 0.78507314
 

The simpson evaluation used the same 50 points with an interval setting of 0.02. As seen from the answer, we can see that it gives a closer answer.

If we use the function 'int' :

Input 6 : int('f',0,1)
Ans. 0.7855690
 

The answer is comparable with simpson method. However, there is some differences between these two methods :

1. The evaluation time for both functions are different. 'int' is much faster than 'simp'.
2. 'simp' has one more setting : an interval setting. Through this setting, we can improve the accuracy of the answer by using a smaller number. Instead, 'int' do not have interval setting.

To improve the accuracy of the integration answer for 'simp' is easy. We can set n = 100 or interval = 0.01 and this is what we got :

Input 7 : simp('f',0,1,0.01)
Ans. 0.78528330
 

It immediately see an improvement to the answer. Actually, as long as we keep on decreasing the interval, the answer would be closer to the exact answer. So, it is considered that this method is a very reliable method. The only disadvantage of this method is that it is slow and need longer evaluation time when the interval get smaller.

For the function 'int', there is also method to improve its accuracy. On idea is to divide 0 to 1 into two partitions and using 'int' to evaluate them separately as follows :

Input 8 : int('f',0,1/2)+int('f',1/2,1)
Ans. 0.78545864
 

We can see that the result is indeed better than just using one partition.

Now, we had seen how the functions of integration worked and we would try to do more examples below :

1. Evaluate ∫0π/2 1/(1 + tan(x)2000) dx 

Input 9 : int('1/(1+pow(tan(x),2000))',0,Pi/2)
Ans. 0.78539816
 

The exact answer for this integral is Pi/4.

2. Evaluate lim n→∞ ∫01 xn/(1 + x2n) dx

To find the limit of an integral, we use function 'table' and evaluate the integral with n set at 1, 10, 100, 1000.

Input 10 : clear('n')
 

Since we had set the value of n previously, we need to make n as a free variable again using clear() command. It not, the table input command 12 that we are going to set will give undefined error. To check if n is a free variable, just type in the symbol n in the input command. If the answer returned is 'n', that means n is a free variable. If a real or complex number returned, then n is not a free variable.

Input 11 : f='pow(x,n)/(1+pow(x,2*n))'
 

Input 12 : table('int(f,0,1)',pow(10,id(4)),'n')
Ans. 0.3466, 0.0703, 0.00787, 4.359e-9
 

The list converges to zero. So, zero is the answer.

3. Evaluate ∫0∞ x e-x2 dx 

Input 13 : f='x*exp(-x*x)'
 

Input 14 : int('f',0,'posinf')
Ans. 0.4993
 

The exact answer for this question is 0.5.

It is also noticed that positive infinity appeared on the upper bound of the integral and it is labelled as 'posinf'. Similarly, for negative infinity, the label is 'neginf'. This integral, where the range or bound of the integration is not finite, is called improper integral.

For the function 'simp', we cannot set infinity as the bound. So when we need to evaluate the above integral using 'simp', we need to find an upper bound for the integral. In this example, we noticed that the value of the function at x = 3, 5 and 10 are :

Input 15 : table('f',(3,5,10))
Ans 0.00037, 6.944e-11, 3.72e-43
 

This means that the function is a very fast decreasing function when x increases. So, we can actually set the upper bound to 3 and it should not change the calculation result too much :

Input 16 : simp('f',0,3,0.06)
Ans. 0.49994
 

In many similar calculation of improper integrals, it is found that 'simp' is actually more accurate than 'int'. Below is an example :

4. Evaluate ∫0∞ 10000 e-0.004x dx

Using 'int' :

Input 17 : f='10000*exp(-0.04*x)'
 

Input 18 : int('f',0,'posinf')
Ans. 179711.7
 

Compare with 'simp' :

Input 19 : simp('f',0,200,0.5)
Ans. 249916.13
 

The exact answer for the question is 250000. So, the answer got from 'int' is not accurate. Although 'simp' is much slower in term of the calculation time than 'int', it is indeed more reliable. Therefore, we suggest two functions are needed when evaluating improper integral. 'int' can be used for the initial fast evaluation and 'simp' for final verification.

Can we improve the accuracy of 'int' by split the integration into two partitions ?

First, we redefine f as a function of 't' and since we do not know where we should split the integration, we define a new function g with the splitting point of the integration as an variable 'x' :

Input 20 : f='10000*exp(-0.04*t)'
 

Input 21 : g='int(f,0,x,t)+int(f,x,posinf,t)'
 

Input 22 : plot('line','g',0,2000,100)
 

Below is the graph showing g as a function of x :

From the graph, we can see that there exists a range of x, (in around 500) that can actually gives the correct answer of 25000 for the integration. So, in this case, splitting can improve the result of 'int' if we choose the right splitting point.

Now we came to double integration.

5. Evaluate ∫02∫01 ex+y dxdy

Actually, the calculation of double integration is not much different from single integration :

Input 23 : f='exp(x+y)'
 

Input 24 : int('int(f,0,1)',0,2,'y')
Ans. 10.978
 

The exact answer for this integral is (e² - 1)(e¹ - 1) :

Input 25 : (exp(2)-1)*(exp(1)-1)
Ans. 10.978199
 

For the above integral, it is x which is being integrated before y but it is equally valid that y can be put to be integrated first and x follows. The equality is shown below :

∫02∫01 ex+y dxdy = ∫01∫02 ex+y dydx

The input commands to evaluate the left handed side formula is shown above with input command 17 and 18. The right handed side formula is entered below and they should give same result :

Input 26 : int('int(f,0,2,y)',0,1)
Ans. 10.978
 

Sometimes, the bounds of a double integral may not be a constant number :

6. Evaluate ∫01∫yy/2 x y dxdy 

Input 27 : f='x*y'
 

Input 28 : int('int(f,y/2,y)',0,1,'y')
Ans. 0.09375
 

The exact value of this integral is 3/32 :

Input 29 : 3/32
Ans. 0.09375
 

7. Evaluate ∫01∫2-x2x2 x1/2 y dydx 

Input 30 : f='sqrt(x)*y'
 

Input 31 : int('int(f,pow(x,2),2-pow(x,2),y)',0,1)
Ans. 0.7621
 

The exact value of this integral is 16/21 :

Input 32 : 16/21
Ans. 0.07619
 

In general, the evaluation time of an integral is proportional exponentially with the level of integration. i.e. suppose 3 seconds is needed to calculate single integration, double integration would takes 9 seconds.

Therefore, one important thing we want to mention is that it is not recommended using 'simp' to evaluate double integral. From our own evaluation, the calculation is too long to get a good answer.

Now, we came to multiple integration. First, we need the latest version of Mathcalc8. At the time of writing, it is version v.2.17. Please install this version or later to continue.

8. Evaluate ∫13∫23∫12 x - y + z dxdydz 

Input 33 : f='x-y+z'
 

Input 34 : g='int(f,1,2)'
 

Input 35 : int('int(g,2,3,y)',1,3,'z')
Ans. 2
 

The exact answer is 2.

9. Evaluate ∫01∫01-x∫01-x-y dzdydx

Input 36 : f='1'
 

Input 37 : g='int(f,0,1-x-y,z)'
 

Input 38 : int('int(g,0,1-x,y)',0,1)
Ans. 0.1667
 

The exact answer is 1/6.

One application of multiple integration is to calculate surface and volume of space that enclosed by curved surface. Below is an example to calculate the 'volume' of a sphere in 4 dimensional space :

We begin with a circle. The area of a circle with radius 1 can be calculated through the following double integration :

10. ∫1-1∫√(1-x2)-√(1-x2) dydx 

Input 39 : f='1'
 

Input 40 : int('int(f,-sqrt(1-x*x),sqrt(1-x*x),y)',-1,1)
Ans. 3.1435342
 

The exact answer is Pi.

Next, we find the volume of a sphere with radius 1 by the following triple integration formula :

11. ∫1-1∫√(1-x2)-√(1-x2)∫√(1-x2-y2)-√(1-x2-y2) dzdydx 

Input 41 : g='int(f,-sqrt(1-x*x-y*y),sqrt(1-x*x-y*y),z)'
 

Input 42 : int('int(g,-sqrt(1-x*x),sqrt(1-x*x),y)',-1,1)
Ans. 4.1913790
 

The exact answer is 4 π /3 :

Input 43 : 4*Pi/3
Ans. 4.1887902
 

As we can see from the two formulas above, their forms looks very similar. It is not difficult to guess a formula for 4D sphere by 'extension' of the form. Below is the formula we guessed :

12. ∫1-1∫√(1-x2)-√(1-x2)∫√(1-x2-y2)-√(1-x2-y2)∫√(1-x2-y2-z2)-√(1-x2-y2-z2) dtdzdydx where

t is the variable name for the 4th dimension. 

Input 44 : clear()
 

This is to ensure that t,z,y,x are free variables. If not, incorrect answer may be resulted.

Input 45 : f='1'
 

Input 46 : g='int(f,-sqrt(1-x*x-y*y-z*z),sqrt(1-x*x-y*y-z*z),t)'
 

Input 47 : h='int(g,-sqrt(1-x*x-y*y),sqrt(1-x*x-y*y),z)'
 

Input 48 : int('int(h,-sqrt(1-x*x),sqrt(1-x*x),y)',-1,1)
Ans. 4.9377010
 

We can compare this value with a formula of the volume of a n dimensional sphere (called (n-1)-sphere) that searched from the Internet :

(2π)n/2 rn / (2 4 ... n) if n is even

2 (2π)(n-1)/2 rn / (1 3 ... n) if n is odd

n being the dimension of the sphere and r is the radius.

According to this formula, the volume of a 4D sphere is (π)²/2 :

Input 49 : Pi*Pi/2
Ans. 4.9348022
 

We can see that our answer matched quite well with the prediction of this formula.