Page 258 - Applied Numerical Methods Using MATLAB

PROBLEMS 247
(d) To investigate how the accident of Jacobian singularity happened, add
h = 1e-5 to the (tentative) solution (rn2) obtained in (c). Does the
result differ from rn2? If not, why? (See Section 1.2.2 and Prob-
lem 1.19.)

>>rn2 + 1e-5 ~= rn2

(e) Charley thought that Jessica just circumvented the Jacobian singularity
problem. To remove the source of singularity, he modified the formula
(5.1.8) into
f((1 + h)x) − f((1 − h)x)

D (x, h) = (P5.5.3)
c2
2hx
and implemented it in another routine “jacob1()” as follows.

function g = jacob1(f,x,h,varargin) %Jacobian of f(x)
if nargin<3, h =.0001; end
h2 = 2*h; N = length(x); I = eye(N);
for n = 1:N
if abs(x(n))<.0001, x(n) =.0001; end
delta = h*x(n);
tmp = I(n,:)*delta;
f1 = feval(f,x + tmp,varargin{:});
f2 = feval(f,x - tmp,varargin{:});
f12 = (f1 - f2)/2/delta; g(:,n) = f12(:);
end

With h = 1e-5 or h = 1e-6 and jacob() replaced by jacob1() in
the routine “newtons()”, type the same statement as in (c) to get a
solution to the problem in Example 4.3 together with its residual error
and check if his scheme works fine.

>>rn3 = newtons(’phys’,1e6,1e-4,100), phys(rn3)

5.4 Numerical Integration of Basic Functions
Compute the following integrals by using the trapezoidal rule, the Simp-
son’s rule, and Romberg method and fill in the following table with the
resulting errors.

2 π/2 1
3
(i) (x − 2x) dx (ii) sin xdx (iii) e −x dx
0 0 0