<< Chapter < Page | Chapter >> Page > |
Sometimes it is rather convenient to use a numerical approach to solve a definite integral. The trapezoid rule allows us to approximate a definite integral using trapezoids.
trapz
CommandZ = trapz(Y) computes an approximation of the integral of Y using the trapezoidal method.
Now, let us see a typical problem.
Given , an analytical solution would produce 39. Use trapz command and solve it
x=2:.1:5;
y=x^2;
. Note the following error prompt:
??? Error using ==>mpower
Inputs must be a scalar and a square matrix.
This is because x is a vector quantity and MATLAB is expecting a scalar input for y. Because of that, we need to compute y as a vector and to do that we will use the dot operator as follows:
y=x.^2;
. This tells MATLAB to create vector y by taking each x value and raising its power to 2.area1=trapz(x,y)
area1 =39.0050
Notice that this numerical value is slightly off. So let us increase the number of increments and calculate the area again:
x=2:.01:5;
y=x.^2;area2=trapz(x,y)
area2 =39.0001
Yet another increase in the number of increments:
x=2:.001:5;
y=x.^2;area3=trapz(x,y)
area3 =39.0000
Determine the value of the following integral:
x=0:pi/100:pi;
y=sin(x);
area1=trapz(x,y)
area1 =1.9998
let us increase the increments as above:
x=0:pi/1000:pi;
y=sin(x);area2=trapz(x,y)
area2 =2.0000
A gas expands according to the law, PV 1.4 =c. Initially, the pressure is 100 kPa when the volume is 1 m 3 . Write a script to compute the work done by the gas in expanding to three times its original volume O. N. Mathematics: 2 by J. Dobinson, Penguin Library of Technology. © 1969, (p. 184) .
Recall that PV diagrams can be used to estimate the net work performed by a thermodynamic cycle, see Wikipedia or we can use definite integral to compute the work done (WD) as follows:
If we rearrange the expression pressure as a function of volume, we get:
By considering the initial state, we can determine the value of c:
From the equation and the equation above, we can write:
By inserting P in WD , we get:
For MATLAB solution, we will consider P as a function of V and WD . Now, let us apply the three-step approach we have used earlier:
v=1:0.001:3;
p=100./v.^1.4;
trapz
function to calculate the work done, the output will be as follows:WorkDone=trapz(v,p)
WorkDone =88.9015
These steps can be combined in an m-file as follows:
clc
disp('A gas expands according to the law, pv^1.4=C')disp('Initial pressure is 100 kPa when the volume is 1 m3')
disp('Compute the work done by the gas in expanding')disp('To three times its original volume')
disp(' ') % Display blank linev=1:.001:3; % Creating a row vector for volume, v
p=100./(v.^1.4); % Computing pressure for volumeWorkDone=trapz(v,p) % Integrating p*dv over 1 to 3 cubic meters
A body moves from rest under the action of a direct force given by where x is the distance in meters from the starting point. Write a script to compute the total work done in moving a distance 10 m. O. N. Mathematics: 2 by J. Dobinson, Penguin Library of Technology. © 1969, (p. 183)
Recall that the general definition of mechanical work is given by the following integral, see Wikipedia :
Therefore we can write:
Applying the steps we followed in the previous examples, we write:
clc
disp('A body moves from rest under the action of a direct force given')disp('by F=15/(x+3) where x is the distance in meters')
disp('From the starting point.')disp('Compute the total work done in moving a distance 10 m.')
disp(' ') % Display blank linex=0:.001:10; % Creating a row vector for distance, x
F=15./(x+3); % Computing Force for xWorkDone=trapz(x,F) % Integrating F*dx over 0 to 10 meters.
The output of the above code is:
A body moves from rest under the action of a direct force given
by F=15/(x+3) where x is the distance in metersFrom the starting point.
Compute the total work done in moving a distance 10 m.WorkDone =21.9951
Z = trapz(Y)
computes an approximation of the integral of Y using the trapezoidal method.Notification Switch
Would you like to follow the 'A brief introduction to engineering computation with matlab' conversation and receive update notifications?