Math 449 - Homework 8 - Solutions
                   
                   
                   

1) Exercise 8 (a, b, c), section 5.5, page 319.

Find the Be'zier curves with degree N and points indicated below.

a) N = 3; {(1,3), (3,-1), (2,4), (3,0)}

P(t) = (1,3)B_03(t) + (3,-1)B_13(t) + (2,4)B_23(t) + (3,0)B_33(t) 

     = (1,3)(1-t)^3 + (9,-3)t(1-t)^2 +(6,12)t^2(1-t) +(3,0)t^3
     
     = ((1-t)^3 + 9t(1-t)^2 + 6t^2(1-t) + 3t^3, 3(1-t)^3 - 3t(1-t)^2 + 12t^2(1-t))
     
     = (1 + 6t - 9t^2 + 5t^3, 3 -12t + 27t^2 - 18t^3)
     
b) N = 4; {(-2, 3), (-1, 3), (3, 5), (3, 4), (2, 3)}

P(t) = (-2, 3)B_04(t) + (-1, 3)B_14(t) + (3, 5)B_24(t) + (3, 4)B_34(t) + (2, 3)B_44(t)

     = (-2, 3)(1-t)^4 + (-4, 12)t(1-t)^3 + (18, 30)t^2(1-t)^2 + 
     
                                           (12, 16)t^3(1-t) + (2, 3)t^4
     
     = (-2 + 4t + 18t^2 - 28t^3 + 10t^4, 3 + 12t^2 - 20t^3 + 8t^4)
     
c) N = 5; {(1, 1), (2, 2), (3, 4), (4, 4), (5, 2), (6, 1)}   

P(t) = (1, 1)B_05(t) + (2, 2)B_15(t) + (3, 4)B_25(t) + (4, 4)B_35(t) + 
                       
                                       (5, 2)B_45(t) + (6, 1)B_55(t) 
                                       
     = (1, 1)(1-t)^5 + (10, 10)t(1-t)^4 + (30, 40)t^2(1-t)^3 + (40, 40)t^3(1-t)^2 + 
                       
                                       (25, 10)t^4(1-t) + (6, 1)t^5     
     
     = (1 + 5t, 1 + 5t + 10t^2 - 30t^3 + 15t^4)          

************************************************************************** 

2) Algorithms and programs 1 and 2, section 5.5, page 319.

The following program produces the Bezier curve for a given set of control

points. The control points are entered as a N+1 by 2 matrix called Points 

whose rows are the coordinate pairs of these points. To use the program

enter 

[X, Y]=bez(Points);

Then plot the curve using 

plot(X, Y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X, Y]=bez(Points)
%Input  - Points is a N+1 by 2 matrix whose rows are the control points
%         of the Bezier curve
%Output - X is the vector of x-coordinates of the Bezier curve
%         Y is the vector of y-coordinates of the Bezier curve

[a b] = size(Points);
N = a-1;
t = 0:0.05:1;
B = (1-t).^N;
X = Points(1,1)*B;
Y = Points(1,2)*B;

for i = 2:N+1
    c = factorial(N)/(factorial(i-1)*factorial(N-i+1));
    B = c*t.^(i-1).*(1-t).^(N-i+1);
    X = X + Points(i,1)*B;
    Y = Y + Points(i,2)*B;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The curves for the control points of problem 8, section 5.5, page 319 are

shown on a separate page.

************************************************************************** 

3) Exercise 4 (a, c), section 6.1, page  335.

Let f(x)=e^x.

(a) We want to compare the approximation to f'(2.3) given by the centered-

formula of order O(h^4) with the value e^(2.3)=9.97418245481472 for h=0.1

and h=0.01.  

Write D_h = [-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)]/12h.

Then D_h = (-exp(2.3 + 2*h) + 8*exp(2.3 + h) - 8*exp(2.3 - h) + exp(2.3 - 2*h))/(12*h)

If h=0.1,  D_0.1  = 9.97414916793670

If h=0.01, D_0.01 = 9.97418245148968

Thus, exp(2.3) - D_0.1 = 3.329 10^-5 and exp(2.3) - D_0.01 = 3.325 10^-9

(c) Using |f^(5)(c)| < exp(2.5), which is approximately 12.18249396, we obtain

the truncation errors: 

E_trunc(f,0.1)  < (0.1)^4*12.18249396/30  = 4.061 10^-5

E_trunc(f,0.01) < (0.01)^4*12.18249396/30 = 4.061 10^-9

************************************************************************** 

4) Exercise 16, section 6.1, page 338.

The function f(x)=ln(x) is given by the following table, which includes

round-off error. The round-off error used was e < 5*10^-6.

       x    |     y
    ________|__________   
    2.8000  |  1.02962
    2.9000  |  1.06471
    2.9980  |  1.09795
    2.9990  |  1.09828
    3.0000  |  1.09861
    3.0010  |  1.09895
    3.0020  |  1.09928
    3.1000  |  1.13140
    3.2000  |  1.16315
    
    
 (a) We  use the formula (-y_2 + 8*y_1 - 8*y_(-1) + y_(-2))/(12*h) to approximate
 
 f'(3.0) for h=0.1 and h=0.001. When h=0.1, 
 
3.0 - 2*h = 2.8       y_(-2) = 1.02962
3.0 -  h  = 2.9       y_(-1) = 1.06471     
3.0 +  h  = 3.1       y_1    = 1.13140
3.0 + 2*h = 3.2       y_2    = 1.16315

so (-y_2 + 8*y_1 - 8*y_(-1) + y_(-2))/(12*h) = 0.33332500000000

For h=0.001, 

3.0 - 2*h = 2.998       y_(-2) = 1.09795
3.0 -  h  = 2.999       y_(-1) = 1.09828    
3.0 +  h  = 3.001       y_1    = 1.09895
3.0 + 2*h = 3.002       y_2    = 1.09928

so (-y_2 + 8*y_1 - 8*y_(-1) + y_(-2))/(12*h) = 0.33583333333342

 (b) We now calculate the total error bounds for each case. 
 
 The total error bound can be estimated by the formula (25, page 328):
 
 |E(f,h)| < 3e/2h + Mh^4/30.
 
 I will use M = 24/(2.8)^5 = 0.14 and e = 5*10^(-6). 
 
 When h=0.1,   we have an error bound 7.5*10^(-5).
 
 When h=0.001, we have an error bound 7.5*10^(-3).

************************************************************************** 

5) Algorithms and Programs 2 (a, b, c), section 6.1, page 338.

The program is written below, with a slight modification so that the 

values of R are also displayed. This is useful since R is an 

approximation for the relative error. I set the tolerance value

at 5*10^(-13). The best I could get is accuracy to 10 significant digits

(using the default parameters).


(a) f(x) = 60*x^45 -32*x^33 +233*x^5 -47*x^2 -77, x=1/sqrt(3)


h              =  0.00000010000000  
derivative     = 75.17349466468202   
relative error =  0.00000000047260


(b) f(x) = tan(cos((sqrt(5) + sin(x))/(1+x.^2))), x=(1+sqrt(5))/3

h              =  0.00000100000000   
derivative     =  1.22859742318626   
relative error =  0.00000000017508


(c) f(x) = sin(cos(1/x)), x=1/sqrt(2)

h              =  0.00000100000000   
derivative     =  1.95155960899807  
relative error =  0.00000000012516

Here is the program:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,n]=difflim(f,x,toler)
%Input   - f is the function input as a string 'f'
%        - x is the differentiation point
%        - toler is the tolerance for the error
%Output  - L=[H' D' E'], where
%        - H is the vector of step sizes
%        - D is the vector of approximate derivatives
%        - E is the vector of error bounds
%        - R is the vector of approximate relative errors
%        - n is the coordinate of the ''best approximation''
max1=15;
h=1;
H(1)=h;
D(1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(1)=0;
R(1)=0;
for n=1:2
    h=h/10;
    H(n+1)=h;
    D(n+1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
    E(n+1)=abs(D(n+1)-D(n));
    R(n+1)=2*E(n+1)/(abs(D(n+1))+abs(D(n))+eps);
end
n=2;
while((E(n)>E(n+1))&(R(n)>toler))&n<max1
    h=h/10;
    H(n+2)=h;
    D(n+2)=(feval(f,x+h)-feval(f,x-h))/(2*h);
    E(n+2)=abs(D(n+2)-D(n+1));
    R(n+2)=2*E(n+2)/(abs(D(n+2))+abs(D(n+1))+eps);
    n=n+1;
end
n=length(D)-1;
L=[H' D' E' R'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

To use the program, I wrote a function file for each function (and called 

the file g.m) then entered the command [L,n]=difflim('g',x,5*10^(-13))

where x is the point where the derivative should be calculated 

(e.g., x=1/sqrt(3) in part (a)).

************************************************************************** 

6) Exercise 4, section 6.2, page 349.

Consider the table for f(x)=cos(x) rounded to four decimal places.

       
       x    |    y
    ________|_________   
    0.9000  |  0.6216
    0.9500  |  0.5817
    1.0000  |  0.5403
    1.0500  |  0.4976
    1.1000  |  0.4536


(a) Find f''(1) using (f_1 - 2*f_0 + f_(-1))/h^2, with h=0.05.

1-h=0.95  f_(-1)=0.5817
1         f_0   =0.5403
1+h=1.05  f_1   =0.4976

The approximate value for f''(1) is  -0.5200 in this case.

(b) Find f''(1) using (f_1 - 2*f_0 + f_(-1))/h^2, with h=0.1.

1-h=0.9   f_(-1)=0.6216
1         f_0   =0.5403
1+h=1.1   f_1   =0.4536

The approximate value for f''(1) is  -0.5400 in this case.

(c) Find f''(1) using (-f_2 + 16*f_1 - 30*f_0 + 16*f_(-1) - f_(-2))/(12*h^2), with h=0.05.

1-2h=0.9   f_(-2)=0.6216
1-h =0.95  f_(-1)=0.5817
1          f_0   =0.5403
1+h =1.05  f_1   =0.4976
1+2h=1.1   f_2   =0.4536

The approximate value for f''(1) is  -0.5133 in this case.

(d) The most accurate answer is (b): f''(1) = -cos(1) = -0.5403.
**************************************************************************