How to solve a definite integral with Scilab

In the article How to calculate a definite integral we have seen how to calculate mathematically a definite integral. Scilab has several built-in functions dedicated to definite integrals calculations. In this tutorial we are going to discuss how to use the function integrate() to calculate definite integrals.

To recall, the definition of a definite integral is:

\[I =\int_{a}^{b} f(x) dx = \left [ F(x) \right ]_a^b= F(b)-F(a)\]

where:

f(x) – continuous function
a – lower limit of integration
b – upper limit of integration
F(x) – antiderivative
I – the result of the integration

\[F(x)=\int f(x) dx\]

In Scilab, the syntax of the integrate() function is:

F=integrate('func', 'x', a, b, AbsTol, RelTol)

where:

'func' – string (Scilab expression) defining the function to be integrated
'x' – string defining the symbolic variable for integration
a – real number (scalar) representing the lower limit of the integration
b – real number (scalar or vector) representing the upper limit of the integration
AbsTol – real number (scalar) representing the absolute error of the numerical integration; default value 1e-8
RelTol – real number (scalar) representing the relative error of the numerical integration; default value 1e-14
F – real number (scalar or vector) representing the result of the integration

The numerical integration is done using a quadrature rule. This is an approximation of the definite integral of a function, as a weighted sum of function values, at specified points within the domain of integration. For example, in the two-point Gauss Quadrature Rule, the integral is approximated as:

\[I =\int_{a}^{b} f(x) dx = \frac{b-a}{2} \cdot \left [ f \left ( \frac{a-b}{2\sqrt{3}} + \frac{b+a}{2} \right ) + f \left ( \frac{b-a}{2\sqrt{3}} +\frac{b+a}{2} \right ) \right ]\]

Let’s take as example the function f(x) defined as:

\[f(x) = 7x^3-8x^2-3x+3\]

and integrate it between a = -0.8 and b = 1.4

Before diving into the calculation process of the integral, let’s do a graphical analysis of the function. We’ll plot the function between the integration limits a and b.

Function f(x) plot with integral areas

Image: Function f(x) plot with integral areas

As you can see, within the integration limits, the function has 3 roots. The result of the integration will be the sum of the areas A, B, C, D. The sum of the positive areas, B and D, is significantly bigger than the sum of the negative areas, A and C. Thus, we expect that the result of the integration to be positive.

Let’s proceed with the integration. First, in order to verify the validity of the integrate function, we’ll solve the symbolic integral, by hand:

\[F(x) = \int {(7x^3-8x^2-3x+3)dx} = \frac{7}{4}x^4-\frac{8}{3}x^3-\frac{3}{2}x^2+3x\]

To find the value of the integral, between the limits a and b, we’ll use a couple of Scilab instructions:

a=-0.8;b=1.4;
deff('y=F(x)','y=(7/4)*x^4-(8/3)*x^3-(3/2)*x^2+3*x');
I = F(b)-F(a);

Now we know the value of the integral I being:

--> I

I =
1.9433333

-->

To get the integration result using the Scilab function integrate(), we’ll use the following Scilab instruction:

I=integrate('7*x^3-8*x^2-3*x+3','x',a,b);

Let’s check the result of the integration:

--> I

I =
1.9433333

-->

As you can see, the result of the Scilab integrate() function is identical with the symbolic calculation.

Alternatively, we can define our f(x) function as a custom Scilab function and use it as an argument for the integrate() function. This approach could be better if the function to be integrated has a complex expression and needs to be defined in a separate file. For our example we are going to use an inline Scilab function definition but the same approach can be applied using *.sci function files.

deff('y=f(x)','y=7*x^3-8*x^2-3*x+3');
I=integrate('f','x',a,b);

Again, checking the result of the integration, gives, as expected, the same value:

--> I

I =
1.9433333

-->

Since working with integrals and Scilab is exciting, let’s integrate the function f(x) in steps, between the limits a, b and the roots, in order to find the individual value for the areas A, B, C and D.

The first step is to find the roots of the function f(x). To do this, we’ll define the function f(x) as a polynomial p(x), using the Scilab function poly():

p=poly([3 -3 -8 7],'x','coeff');

--> p

p =
         2   3
3 -3x -8x +7x

-->

To find the roots of the polynomial, we need to call the roots() function with the argument p:

r=roots(p);

--> r

r =
-0.6276878
1.2029662
0.5675787

-->

The vector r contains the roots of our function but as complex variables. Let’s convert them to real variables, and sort them from minimum to maximum, using the Scilab real() and gsort() functions:

r=gsort(real(r),'g','i');

--> r

r =
-0.6276878
0.5675787
1.2029662

-->

Next, we’ll create a vector iLim which starts from the a limit, contains the roots r and stops at the b integration limit:

iLim = [a r' b];

--> iLim

iLim =
-0.8 -0.6276878 0.5675787 1.2029662 1.4

-->

With a Scilab FOR loop we can integrate our function f(x) in order to calculate the areas A, B, C and D individually:

for i=1:length(iLim)-1
   X(i)=integrate('f','x',iLim(i),iLim(i+1));
end

--> X

X =
-0.2650555
2.4564698
-0.4527833
0.2047023

-->

Each value of the vector X represents the area of each region A, B, C and D (in this order). As expected A and C are negative, B and D are positive and B has the highest value. These results match with the graphical representation of the f(x) depicted in the image above.

If we sum up all the values of the vector X we should get the same value of I = 1.9433333.

--> sum(X)

ans =
1.9433333

-->

The Scilab script below contains all the Scilab instructions mentioned in the tutorial together with the instruction for generating the plot of the function f(x).

clear, clc
// Define the function f(x) and integration limits
deff('y=f(x)','y=7*x^3-8*x^2-3*x+3');
a=-0.8;b=1.4;

// Define the antiderivative F(x) and evaluate the integral
deff('y=F(x)','y=(7/4)*x^4-(8/3)*x^3-(3/2)*x^2+3*x');
I=F(b)-F(a);

// Integrate the expression of the function f(x)
I=integrate('7*x^3-8*x^2-3*x+3','x',a,b);

// Integrate the custom function f(x)
I=integrate('f','x',a,b);

// Define the polynomial p(x), find and sort the roots
p=poly([3 -3 -8 7],'x','coeff');
r=roots(p);
r=gsort(real(r),'g','i');

// Calculate the integral for each area: A, B, C and D
iLim = [a r' b];
for i=1:length(iLim)-1
 X(i)=integrate('f','x',iLim(i),iLim(i+1));
end
I = sum(X);

// Plot the f(x) function, roots and intgeration areas
x=a:0.01:b;
plot(x,f(x),'LineWidth',2), xgrid();
xlabel('x','FontSize',3);
ylabel('$f(x)=7x^3-8x^2-3x+3$','FontSize',4);
title('x-engineer.org','FontSize',2,'Color','blue')
xfpoly(x,f(x));
h1=get("hdl"); 
h1.data=[x' f(x)']
h1.mark_mode="off";
h1.polyline_style=6;
h1.bar_width=0.007;
h1.foreground=4;
h1.background=4;
plot(x,f(x),'LineWidth',2) // overlap grid on poligon
plot(r,[0 0 0]')
h2=get("hdl");
h2.children.line_mode='off';
h2.children.mark_mode='on';
h2.children.mark_style=9;
h2.children.mark_size=8;
h2.children.mark_foreground=10;
h2.children.thickness=2;
xstring(-0.77,-0.8,'A')
xstring(-0.15,1.4,'B')
xstring(0.89,-0.7,'C')
xstring(1.3,0.3,'D')

For any questions, observations and queries regarding this article, use the comment form below.

Don’t forget to Like, Share and Subscribe!

Leave a Reply

Ad Blocker Detected

Dear user, Our website provides free and high quality content by displaying ads to our visitors. Please support us by disabling your Ad blocker for our site. Thank you!

Refresh