Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Numerical integration usibg Matlab

Subject: Numerical integration usibg Matlab

From: Jules Kassmann

Date: 12 Mar, 2011 12:13:05

Message: 1 of 6

Hello,

I am trying to calculate the power produced by a gas turbine.
It is equal to the opposite of the integration of V dP between the inlet and the outlet.
As it is about a gas (assumed ideal), we know that V = (nRT)/P
To calculate the power, I have to use the build-in 'quad' function, and to implement the trapezoidal rule.
For the trapezoidal rule, I first need to let the user choose the number of sub-intervals, and then calculate the percentage error with the quad method.
I also need to plot the percentage error if [10, 20, 30, 40, 60, 80, 100] sub-intervals are used.

Here is my code:

_______________________________________________________________________

n = 0.2; %Molar flow rate, in kmol/sec
R = 8.314; %Universal gas constant, in kj/kmol K
T = 500; %Temperature, in K
Pi = 600; %Inlet pressure, in kPa
Po = 150; %Outlet pressure, in kPa


W1 = quad(@(P) (-n*R*T)./P ,Pi,Po) %Power found using the quad function


t = input('Choose the number of pressure points (more points => more accurate integration: ');
                                                           %Let the user choose the number of pressure points for the integration
x1 = linspace(Pi,Po,t);
y1 = 1./x1; %Function of the pressure
W2 = sum(diff(x1).*(y1(1:length(y1)-1)+diff(y1)/2)) %Power found using the trapezoidal rule


difference1 = (abs(W1-W2)/W1)*100; %Percentage difference between the value foud using the quad function and the trapezoidal rule
fprintf('There is a %3.2f percent difference between the quadrature and the trapezoid reults',difference1)



N = [10,20,30,40,60,80,100]; %Array containing the number of pressure points for the integration

for n = N(1:length(N))
    x2 = linspace(Pi,Po,n);
    y2 = 1./x2; %Function of the pressure
    W(n) = sum(diff(x2).*(y2(1:length(y2)-1)+diff(y2)/2)); %Power found using the trapezoidal rule for each number of pressure points
    difference(n) = (abs(W1-W(n))/W1)*100 %Percentage difference between the value foud using the quad function and the trapezoidal rule
end

plot(difference,N)

_______________________________________________________________________


There is obviously an error in the W(n) function.
Thank you for your help.

Jules

Subject: Numerical integration usibg Matlab

From: Roger Stafford

Date: 12 Mar, 2011 17:30:21

Message: 2 of 6

"Jules Kassmann" <jules.kassmann@gmail.com> wrote in message <ilfo0h$1tl$1@fred.mathworks.com>...
> ........
> y2 = 1./x2;
> .......
- - - - - - - -
  When using 'quad' you have multiplied by -n*R*T but not when doing the trapezoidal integration, where you just took the reciprocal.

Roger Stafford

Subject: Numerical integration usibg Matlab

From: Jules Kassmann

Date: 12 Mar, 2011 18:11:04

Message: 3 of 6

"Roger Stafford" wrote in message <ilgajd$dqt$1@fred.mathworks.com>...
> "Jules Kassmann" <jules.kassmann@gmail.com> wrote in message <ilfo0h$1tl$1@fred.mathworks.com>...
> > ........
> > y2 = 1./x2;
> > .......
> - - - - - - - -
> When using 'quad' you have multiplied by -n*R*T but not when doing the trapezoidal integration, where you just took the reciprocal.
>
> Roger Stafford


Thank you.
That solves the problem of the first calculation. After that, in the for loop, the remaining problem is that I get 10, then 20, then 30, then 40, then 60, then 80 and finally 100 values (so overall 340 values). most of them are zeros...
I understand that I get the number of values of my N array [10, 20, 30, 40, 60, 80, 100] but I have no idea why.

Subject: Numerical integration usibg Matlab

From: Jules Kassmann

Date: 12 Mar, 2011 19:09:03

Message: 4 of 6

I have spotted a beginner mistake: n was reallocated during the for loop.
Now I have this piece of code containing an error:

________________________________________________
n = 0.2; %Molar flow rate, in kmol/sec
R = 8.314; %Universal gas constant, in kj/kmol K
T = 500; %Temperature, in K
Pi = 600; %Inlet pressure, in kPa
Po = 150; %Outlet pressure, in kPa


W1 = quad(@(P) (-n*R*T)./P ,Pi,Po) %Power found using the quad function

N = [10,20,30,40,60,80,100]; %various values of sub-intervals for the trapezoid rule

for t = N(1:length(N)) %for each value of N,
    x2 = linspace(Pi,Po,t); %I define an array for the integration
    y2 = 1./x2;
    W(t) = (-n*R*T)*sum(diff(x2).*(y2(1:length(y2)-1)+diff(y2)/2))
    difference(t) = (abs(W1-W(t))/W1)*100
                                               %percentage difference between quad and trapezoid
end

plot(difference,N)
__________________________________________________
I think the error is in the expression of W(t).
I also don't know if I have to put x2(t) or x2 alone (same with the other variables of the for loop.

Subject: Numerical integration usibg Matlab

From: Roger Stafford

Date: 12 Mar, 2011 19:51:04

Message: 5 of 6

"Jules Kassmann" <jules.kassmann@gmail.com> wrote in message <ilggcf$n0a$1@fred.mathworks.com>...
> ......
> for t = N(1:length(N)) %for each value of N,
> x2 = linspace(Pi,Po,t); %I define an array for the integration
> y2 = 1./x2;
> W(t) = (-n*R*T)*sum(diff(x2).*(y2(1:length(y2)-1)+diff(y2)/2))
> difference(t) = (abs(W1-W(t))/W1)*100
> end
> ......
- - - - - - - - -
  The error lies in your use of 't'. To use it as an index into 'W' and 'difference' it should consist of successive integers, not values of N. You should write:

for t = 1:length(N)
 x2 = linspace(Pi,Po,N(t)); % <-- Use N(t) here
 y2 = ...
 W(t) = ...
 difference(t) = ...
end

  One additional remark: I don't understand why you used 'quad' as a means of comparison with your trapezoidal results when the integral of 1/P has a known analytic solution, namely the natural logarithm of P. Why not compare your W results with -n*R*T*log(Po/Pi)?

Roger Stafford

Subject: Numerical integration usibg Matlab

From: Jules Kassmann

Date: 12 Mar, 2011 20:54:04

Message: 6 of 6

"Roger Stafford" wrote in message <ilgir8$b5$1@fred.mathworks.com>...
> "Jules Kassmann" <jules.kassmann@gmail.com> wrote in message <ilggcf$n0a$1@fred.mathworks.com>...
> > ......
> > for t = N(1:length(N)) %for each value of N,
> > x2 = linspace(Pi,Po,t); %I define an array for the integration
> > y2 = 1./x2;
> > W(t) = (-n*R*T)*sum(diff(x2).*(y2(1:length(y2)-1)+diff(y2)/2))
> > difference(t) = (abs(W1-W(t))/W1)*100
> > end
> > ......
> - - - - - - - - -
> The error lies in your use of 't'. To use it as an index into 'W' and 'difference' it should consist of successive integers, not values of N. You should write:
>
> for t = 1:length(N)
> x2 = linspace(Pi,Po,N(t)); % <-- Use N(t) here
> y2 = ...
> W(t) = ...
> difference(t) = ...
> end
>
> One additional remark: I don't understand why you used 'quad' as a means of comparison with your trapezoidal results when the integral of 1/P has a known analytic solution, namely the natural logarithm of P. Why not compare your W results with -n*R*T*log(Po/Pi)?
>
> Roger Stafford


Thank you very much, now everything works smoothly =)

I know the example is quite simple, but it is a textbook exercise and the question was to use the quad function. (MATLAB for Engineers by Holly Moore, Problem 12.19 modified for the purpose of the lecture)


Jules

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us