Page 79 - Applied Numerical Methods Using MATLAB

68 MATLAB USAGE AND COMPUTATIONAL ERRORS

function INTf = trpzds(f,a,b,N)
%integral of f(x) over [a,b] by trapezoidal rule with N segments
if abs(b - a) < eps | N <= 0, INTf = 0; return; end
h = (b - a)/N; x = a+[0:N]*h;
fx = feval(f,x); %values of f for all nodes
INTf = h*((fx(1)+ fx(N + 1))/2 + sum(fx(2:N))); %Eq.(5.6.1)

1.24 Parameter Passing through varargin

Consider the integration routine ‘trpzds(f,a,b,N)’ in Section 5.6. Can
you apply it to compute the integral of a function with some parame-
ter(s), like the ‘Bessel_integrand(x,beta,k)’ that you defined in Prob-
lem 1.21? If not, modify it so that it works for a function with some param-
eter(s) (see Section 1.3.6) and save it in the M-file named ‘trpzds_par.m’.
Then replace the ‘quad()’ statement in the program ‘nm1p21b’ (introduced
in P1.21) by an appropriate ‘trpzds_par()’ statement (with N = 1000)and
run the program. What is the discrepancy between the integration results
obtained by this routine and the recursive computation based on Problem
1.21.4(b)? Is it comparable with that obtained with ‘quad()’? How do you
compare the running time of this routine with that of ‘quad()’? Why do
you think it takes so much time to execute the ‘quad()’ routine?
1.25 Adaptive Input Argument to Avoid Runtime Error in the Case of Different
Input Arguments
Consider the integration routine ‘trpzds(f,a,b,N)’ in Section 5.6. If some
user tries to use this routine with the following statement, will it
work?
trpzds(f,[a b],N) or trpzds(f,[a b])

If not, modify it so that it works for such a usage (with a bound vector as
the second input argument) as well as for the standard usage and save it in
the M-file named ‘trpzds_bnd.m’. Then try it to find the intergal of e −t
for [0,100] by typing the following statements in the MATLAB command
window. What did you get?
>>ftn=inline(’exp(-t)’,’t’);
>>trpzds_bnd(ftn,[0 100],1000)
>>trpzds_bnd(ftn,[0 100])
1.26 CtFT(Continuous-Time Fourier Transform) of an Arbitrary Signal
Consider the following definitions of CtFT and ICtFT(Inverse CtFT) [W-4]:

∞
X(ω) = F{x(t)}= x(t)e −jωt dt: CtFT (P1.26.1a)
−∞
1 ∞ jωt
−1
x(t) = F {X(ω)}= X(ω)e dω: ICtFT (P1.26.1b)
2π
−∞