(8.11) The files for this problem are contained in the folder s12c7p2 matlab.
The Matlab script is:
1function s12c7p2
2clc
3close all
14 I=2*(f(1)+f(2))/2;
15 fprintf(‘Trapezoidal rule: \t I = %8.6f\n’,I)
16
17 %simpsons 1/3 rule
18 f = get points(3);
29 I=2*(7*f(1)+32*f(2)+12*f(3)+32*f(4)+7*f(5))/90;
30 fprintf(‘Booles rule: \t\t I = %8.6f\n’,I)
31
32 %5th order
33 f = get points(6);
44 I = ntrapezoid(ntraps);
45 fprintf(‘100 trapezoidal rules: \t I = %8.6f\n’,I)
46
47 %Gauss quadrature
48 I = feval(-1/sqrt(3))+ feval(1/sqrt(3));
58 I two = ntrapezoid(100);
59 I = I two + (I two–I one)/99;
60 fprintf(‘Richardson (10/100): \t I = %8.6f\n’,I)
61
62
71 out = 2*fsum/(2*ntraps);
72
73
74 function out = get points(n)
75 %does the function evaluations for n evenly spaced points in …
85 out = exp(-3*x)*log(2+sin(pi*x));
The output of this program is
Trapezoidal rule: I = 13.956743
Simpsons 1/3 rule: I = 5.576444
Simpsons 3/8 rule: I = 4.036023
(8.12) The files for this problem are contained in the folder awc3 matlab
The Matlab script is:
490 Interpolation and integration
1function awc3
2clc
12 U = zeros(n);
13
14 j = 1;
15 for i = 1:n
16 U(i,1) = (f(i+1) – f(i))/(x(i+j) – x(i));
27 coeff(1,2:n+1) = U(1,:);
28 end
29
30 function y = integrate interpolate f(x,xpts,coeff,n) %find y …
value of interpolation; program 8.2
41 function out = integrate multiple trapezoidal(xmin,xmax,nt) …
%mutliple trapezoid integration; program 8.5
42 x = linspace(xmin,xmax,nt+1);
43 I = getf(xmin) + getf(xmax);
44 if nt >1
54 [n,coeff] = integrate interpolation tabulated(f,r); %find the …
coefficients
55
56 a = 0:0.01:0.8; %create x axis for better resolution for plotting
57 y int = zeros(1,81); %matrix to store interpolated polynomial …
65 legend(‘Original’,‘Interpolated’,‘Location’,‘NorthEast’)
66 saveas(g,‘awc3 solution fig.eps’,‘psc2’)
67
68 %%%%%
69
79 fprintf(‘\t b %1d = %10.6f \t x %1d = %10.6f …
\n’,m-1,coeff(m),m-1,f(m))
80 end
81 display(tau)
82
0.5
1
1.5
2.5
3
3.5
4
Original
Interpolated
492 Interpolation and integration
The text output of the program is
1Interpolation coefficients and points are:
14 3.1215
(8.13) The files for this problem are contained in the folder awc2 matlab
The Matlab script is:
1function awc2
2clc
3close all
13 j = 1;
14 for i = 1:n
15 U(i,1) = (f(i+1) – f(i))/(x(i+j) – x(i));
16 end
17
28
29 function y = integrate interpolate f(x,xpts,coeff,n) %find y …
value of interpolation; program 8.2
30 y = coeff(1);
31 for i = 2:n+1
integrate multiple trapezoidal(xmin,xmax,nt,side) …
%mutliple trapezoid integration; program 8.5
41 x = linspace(xmin,xmax,nt+1);
42 if side == 1
43 I = getf top(xmin) + getf top(xmax);
54 end
55 end
56 else
57 display(‘1 or 2 only for side’)
58 end
68 y int = zeros(1,91); %matrix to store interpolated polynomial …
points
69 for b = 1:91
70 y int(b) = integrate interpolate f(a(b),t,coeff,n); …
%evaluate polynomial at different points
80 %%%%%
81 function out = getf top(x) %match top integrand
82 out = x *integrate interpolate f(x,t,coeff,n);
494 Interpolation and integration
83 end
93 fprintf(‘Interpolation coefficients and points are:\n’)
94 for m = 1:n+1
95 fprintf(‘\t b %1d = %10.6f \t x %1d = %10.6f …
\n’,m-1,coeff(m),m-1,t(m))
96 end
t (s)
0 100 200 300 400 500 600 700 800 900
c (ppm)
-5
0
5
10
15
20
25
Original
Interpolated
The text output of the program is
(8.14) The files for this problem are contained in the folder awc1 matlab
The Matlab script is:
1function awc1
2clc
3close all
13
14 j = 1;
15 for i = 1:n
16 U(i,1) = (f(i+1) – f(i))/(x(i+j) – x(i));
17 end
28 end
29
30 function y = integrate interpolate f(x,xpts,coeff,n) %find y …
value of interpolation; program 8.2
31 y = coeff(1);
%mutliple trapezoid integration; program 8.5
42 x = linspace(xmin,xmax,nt+1);
43 I = getf(xmin) + getf(xmax);
44 if nt >1
45 for i = 2:nt
coefficients
55
56 a = 12:.05:19; %create x axis for better resolution for plotting
57 y int = zeros(1,141); %matrix to store interpolated polynomial …
points
66 legend(‘Original’,‘Interpolated’,‘Location’,‘SouthEast’)
67 saveas(g,‘awc1 solution fig.eps’,‘psc2’)
68 %%%%%
69
70 function out = getf(x)
79 for m = 1:n+1
80 fprintf(‘\t b %1d = %10.6f \t x %1d = %10.6f …
\n’,m-1,coeff(m),m-1,t(m))
The output file is:
t (hr)
12 13 14 15 16 17 18 19
c (ppm)
-500
0
500
1000
1500
2000
2500
3000
Original
Interpolated
The text output of the program is
1Interpolation coefficients and points are:
20 0.8295
(8.15) The files for this problem are contained in the folder s14c7p3 matlab The plot os-
cillates because the integral is approximated by a number of trapezoids, which alter–
nately leads to larger and smaller areas as you increase the number of trapezoids. The
amplitude of the oscillations are decreasing and they converge to a constant value.
Functions having sharp peaks are in general troublesome. Using trapezoidal rule,
498 Interpolation and integration
0 10 20 30 40 50 60
−6
−4
−2
0
2
4
6
Number of intervals
Amount adsorbed
1function out = hw7p3()
10 z i = 0.5;
11 MaxNseg = 60;
12 vIntegral = zeros(MaxNseg,1);
13
14 for i = 1:MaxNseg
25 end
26 %End of main program
27
28
29 %Subfunction to implement multiple trapezoidal rule.
38
39
40 %Subfunction to evaluate the integrand.
41 %Input: vector X
42 %Output: vector Y, where Y(i) = C(z(i)) – 1
3.32, we already showed that the liquid volume is 0.3529 and the vapor volume is
1.5221. You also need the intermediate volume. I used the Newton’s method code
from Problem 3.32 with an initial guess of 0.7 to get the intermediate volume of
0.6626. The two integrals that should be equal are
The two integrals are close but not exactly equal. The Newton’s method scheme we
asked you to use on Problem 3.32 is much more accurate because you solve the
original non-linear equations. In the integral construct, you have to deal with the
error in the numerical integration.
The Matlab script is:
500 Interpolation and integration
1function s12c7p3
2clc
3close all
13
14
15 %get each of the roots to set the bounds for the integrals
16 Vl = find V(P,0.3,T,a,b);
17 Vmid = find V(P,0.7,T,a,b);
%6.4f\n’,Vl,Vmid)
32 x = gauss quad convert(Vl,Vmid);
%6.4f\n’,Vmid,Vg)
38 x = gauss quad convert(Vmid,Vg);
39 Ivdw = (Vg-Vmid)/2*(vdw(x(1),T,a,b)+vdw(x(2),T,a,b)); …
%integral of vdw P
40 I = Ivdw – P*(Vg-Vmid); %compute the area under the curve
Interpolation and integration 501
51 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
52
53 function out=plot vdw(T,a,b,Vmin,Vmax,Pmin,Pmax)
54 %Plot vdw EOS between Vmin,Vmax and Pmin,Pmax. returns h to …
64 ylabel(‘$P$ (MPa)’,‘FontSize’,14)
65 top = strcat(‘vdw Equation for $T = $’,num2str(T),‘K’);
66 title(top,‘FontSize’,14)
67 out = h;
68
77 function out = find V(P,V,T,a,b)
78 %executes Newton’s method to find a value of V in L/mol when …
the vdw
79 %equation of state is written as f(V) = 0 with the pressure in …
MPa, the
90 fprintf(‘Failed to compute volume for Vinit = %6.4f …
‘,Vinit)
91 fprintf(‘P = %6.4f and T = %6.2f’,P,T)
92 fVeval = 0; %stop the loop
93 V = 0; %write nonsense for the answer
102 out = R*T/(V-b) – a/Vˆ2 – P;
103
104 function out = dfvdw(V,T,a,b)
105 %computes the value fo the derivative f’ of the vdw equation …
in the form
Interpolation and integration 503
For keven, we have
n−1
X
2
i=0
If we order the unknowns as z=[ci;xi], then the entries to the Jacobian entries for
k=0 : n−1 are the derivatives with respect to c,
For k=n: 2n−1, the derivatives are with respect to x
See the Matlab program output for the coefficients for 5th order polynomials. You
need a very good guess (which helps if you looked up the answer somewhere, since
this result is known) and you need to be careful not to divide by zero. Also, since
Matlab counts from 1, not from 0, you need to be careful about the even/odd values
of k.
1function s13c7p12
2clc
3close all
4set(0,‘defaulttextinterpreter’,‘latex’)
13 %set up the initial guess (this is a good guess)
14 c = [0.7;0.5;0.5;0.5;0.7];
15 x = [-0.9;-0.5;0;0.5;0.9];
16
17 %compile into a single vector
504 Interpolation and integration
28 z = z + delta;
29 R = getR(z,n);
30 k = k+1;
31 err = norm(R);
42 fid = fopen(‘homework07p2 text.txt’,‘w+’); %for file
43 fprintf(‘\n\nThe coefficients and evaluation points are\n’)
44 fprintf(fid,‘The coefficients and evaluation points are\n’); …
%for file
45 for i = 1:n
%8.6f\n’,i-1,c(i),x(i));
47 fprintf(fid,‘i = %3d \t c i = %8.6f \t x i = …
%8.6f\n’,i-1,c(i),x(i)); %for file
48 end
49 fclose(fid);
50
51
62 for j = 1:n
63 R(i,1) = R(i,1) + c(j)*x(j)ˆk;
64 end
65 end
66 %write the constant parts
76 c = z(1:n);
77 x = z(n+1:2*n);
78
79 J = zeros(2*n,2*n);
Interpolation and integration 505
90 J(i,m) = 0; %need to avoid NAN with x = 0
91 end
92 end
93 end
94
95 out = J;
The output to the screen is
1k = 0 err = 1.426358e+00
2k = 1 err = 9.518770e-03