4.12

Evaluate and interpret the condition numbers for

(a) f(x) = sqrt(abs(x-1)) +1       at x=1.00001
(b) f(x) = e ^ -x                          at x=10
(c) f(x) = sqrt(x^2 +1) - x          at x=300
(d) f(x) = (e^-x - 1)/x                 at x=0.001
(e) f(x) = sin(x)/1+cos(x)         at x=1.0001pi


Matlab Code:

e = exp(1);

%(a) f(x) = sqrt(abs(x-1)) +1

xa = -2:.001:2;

fa = sqrt(abs(xa-1)) +1;

derfa = (-1+xa)/(2 .*((-1+xa).^2).^(3/4));

cna = xa .* derfa ./fa;

figure(1);

hold on

plot(xa,fa)

plot(xa,cna, 'r')

hold off

%(b) f(x) = e ^ -x

xb = -10:1:10;

fb = e .^ -xb;

derfb = -e .^ -xb;

cnb = xb .* derfb ./fb;

figure(2);

hold on

plot(xb,fb)

plot(xb,cnb, 'r')

hold off

%(c) f(x) = sqrt(x^2 +1) - x

xc = 0:10:300;

fc = sqrt(xc .^2 +1) - xc;

derfc = -1+xc ./sqrt(1+xc .^2);

cnc = xc .* derfc ./ fc;

figure(3);

hold on

plot(xc, fc)

plot(xc,cnc,'r')

hold off

%(d) f(x) = (e^-x - 1)/x;

xd = -.5:0.001:.05;

fd = (e.^-xd - 1)./xd;

derfd = -e .^(-1 - xd);

cnd = xd .* derfd ./ fd;

figure(4);

hold on

plot(xd, cnd, 'r')

plot(xd,fd)

hold off

%(e) f(x) = sin(x)/1+cos(x);

xe = -1 * pi: .0001 * pi: 2 * pi;

fe = sin(xe) ./(1+cos(xe));

derfe = 1 ./(1+cos(xe));

cne = xe .* derfe ./ fe;

figure(5);

hold on

plot(xe,fe)

plot(xe, cne, 'r')

hold off

The results of the code are these graphs:

A:

B:

C:

D:

E:

As shown, the cn(conditional number) for a is about -.25, meaning the error is attenuated. The cn for b is about 0, so it is attenuated. The cn for c is -1, meaning the relative error is identical in x and f(x). The cn for d is about 0, so it is attenuated. The cn for e appears to be 0 at first, but upon closer inspection, 

it becomes apparent that at ~3.1419, the point at which we are evaluating the condition number, there is a discontinuity which goes to infinity, showing that at that point, the error is magnified enormously

4.8

The Stefan-Boltzmann law can be employed to estimate the rate of radiation of energy H from a surface, as in
 H =  AeoT^4
where H is in watts, A = the surface area (m^2), e = the emissivity that characterizes the emitting properties of the surface (dimensionless), o = a universal constant called the Stefan-Boltzmann constant (=5.67 x 10^-8 W m^-2 K^-4), and T= absolute temperature (k). Determine the error of H for a steel plate with A = .15 m^2, e = .90, and T = 650 + or - 20. Compare your results with the exact error. Repeat the computation but with T = 650 + or - 40.  Interpret your results.

o = 5.67 *10^-8;

A = 0.15;

e = 0.9;

T = 650 + (-20:20);

%T = 650 + (-40:40);

H = @(T) A* e * o * T.^4;

hold on

plot(T, H(T))

hold off

This code yields the graph of H showing a range of error of approximately 350 degrees with 

T = 650 ± 20

As would stand to reason, when we double the range of uncertainty for T to 40(T=650±40), the range of error about doubles

4.5

 Use zero- through third-order Taylor series expansions to predict f(3)
 for

 f(x) = 25*x^3 - 6*x^2 +7*x - 88

using a base point at x = 1

x = 0:0.01: 3.5;

f_x = @(x) 25*x.^3 - 6 * x.^2+7*x -88;

f1_x = @(x) 3* 25*x.^2 - 2* 6 * x.^1+7;

f2_x = @(x) 2* 3* 25*x - 2* 6;

f3_x = @(x) 2* 3* 25;

T_0_1 = f_x(1) * x.^0;

T_1_1 = T_0_1 + f1_x(1) * (x-1).^1;

T_2_1 = T_1_1 + f2_x(1) .* (x-1).^2/2;

T_3_1 = T_2_1 + f3_x(1) * (x-1).^3/6;

hold on

plot(x,f_x(x))             %blue

plot(x, T_0_1, 'r')       %red

plot(x, T_1_1, 'm')     %magenta

plot(x, T_2_1, 'y')      %yellow

%plot(x, T_3_1, 'g')   %green

hold off

With the third order Taylor plot hidden, the graph looks like 

Showing that with increasing orders of the taylor series, the curve more closely approximates the original equation, until we show the third order equation, at which point

one can see that the third order taylor exactly follows the original equation. This stands to reason, as it is mentioned in the text that in general, the nth order Taylor series will be exact for an nth order polynomial, and the equation we are approximating has a largest order term of 25x^3, making it 3rd order. 

4.3

Use the Taylor series to estimate f(x) = e^-x


Matlab code:

a = 0.2;

x = -5: .01 : 5;

f =@(a) exp(-a);

f_1 =@(a) -exp(-a);

f_2 =@(a) exp(-a);

f_3 =@(a) -exp(-a);

T_0 = f(a) * x.^0;

T_1 = T_0 + f_1(a) * (x - a);

T_2 = T_1 + ((f_2(a))/factorial(2)) *(x - a) .^2;

T_3 = T_2 + ((f_3(a))/factorial(3)) *(x - a) .^3;

hold on

plot(x,exp(-x))

plot(x,T_0, 'r')

plot(x,T_1,'g') 

plot(x,T_2, 'k')

plot(x,T_3, 'm')

hold off

This plots the expansions out to the third order as requested by the problem and produces the graph:

As you can see, the zeroth order taylor in red is the farthest from approximating the actual equation in blue, the first order in green is closer, the second in black even closer, and the third in magenta is closest. They all converge at .2, the x we were solving for. 

3.11

Use 5-Digit arithmetic with chopping to determine the roots of the following Eqs (3.12) and (3.13)


3.12: 
x = (-b± sqrt(b^2-4ac))/2a


3.13
x=-2c/(b± sqrt(b^2-4ac))


x^2 - 5000.002x+10


Compute percent relative errors for your results

a = chop(1,5);

b = chop(-5000.002,5);

c = chop(10,5);

  %3.12

  

  x1 = chop((-b + sqrt(b^2 - 4*a*c))/(2*a),5);

  x2 = chop((-b - sqrt(b^2 - 4*a*c))/(2*a),5);

  

  

  %3.13

  

  x1b = chop((-2*c)/(b + sqrt(b^2 - 4*a*c)),5) ;  

  x2b = chop((-2*c)/(b - sqrt(b^2 - 4*a*c)),5) ; 

  

  display(x1);

  display(x1b);

  display(x2);

  display(x2b);

  

  percenterror = ((5000 - x1)/5000) *100;

  percenterrorb = ((5000 - x1b)/5000) *100;

  percenterror2 = ((.002 - x2)/.002) *100;

  percenterror2b = ((.002 - x2b)/.002) *100;

The roots are 5000 and .002

3.6

Evaluate e^-5 using two approaches

 Equation 1
e^ -x = 1 - x + x^2/2 - x^3/3! +......
 Equation 2
1/(e^ -x = 1 + x + x^2/2 + x^3/3! +......)

The actual value of e^5 is 6.737947 * 10^-3

Compute error of each method for 20 terms and compare

n = (0:19); 

x= 5; 

signs = (-1).^n; 

x2n = x.^n ; 

enx = cumsum(signs.* x2n./factorial(n));% eq1

epx = 1../(cumsum(x.^n ./factorial(n)));% eq2

m=1:.001:10;

true = 6.737947 * 10^-3 ;

hold on

plot(enx); %eq1

plot (epx, 'r');% eq2(red)

plot(m,true); 

hold off

The equation in blue is the first equation, the second is in red. They both converge around 15, but the second gets close sooner

3.5

The infinite series  converges on a value of f(n) = pi^4/90 as n approaches infinity. Compute the equation from i = 1 to 10,000 and 10,000 to 1. 

Pi^4/90 calculates to 1.08

syms k;

s1 = symsum( 1/k^4,1,10000);

ezplot(s1)

I plotted it again using

s2 = symsum( 1/k^4,10000,1);

 and the graph was a line on 0, I'm not sure why it doesn't converge the same way

3.2

%71563 from base 8 to base 10

x = base2dec('71563',8)

%3.14 from base 8 to base 10

y = 3 * 8^0 + 1 * 8^-1 + 4 * 8^-2

x =

       29555

y =

    3.1875

2.14

Two distances are required to specify the location of a point relative to an origin in two dimensional space:

• The horizontal and vertical distances (x,y) in Cartesian coordinates
• The radius and angle (r, theta) in radial coordinates

It is relatively straightforward to compute Cartesian coordinates (x,y) on the basis of polar coordinates (r,theta). The reverse process is not so simple. The radius can be computed by the following formula: 
r = sqrt(x^2 + y^2)
If the coordinates lie within the first and fourth coordinates (i.e., x>0), then a simple formula can be used to compute theta:
theta = arctan(y/x)
The difficulty arises for the other cases. The following table summarizes the possibilities:

Write a program to calculate r and theta and test is by filling in this table:

x= [1,1,0,-1,-1,-1,0,1,0];

y= [0,1,1,1,0,-1,-1,-1,0];

theta = atan2(y,x);

theta = theta +(theta<0)*2*pi

r = (sqrt(x.^2+y.^2))

2.9

Economic formulas are available to compute annual payments for loans. Suppose that you borrow an amount of money P and agree to repay it in n annual payments at an interest rate of i. The formula to compute the annual payment A is A = P * ((i (1+i)^n)/((1+i)^n - 1)). Write a program to compute A. Test it with P = $55000 an an interest rate of 66%. Compute results for n = 1,2,3,4,5.


Matlab Code:

p = 55000;

i = .066;

n = 1:5;

A = p * ((i*(1+i).^n)./((1+i).^n - 1));

Results:

2.8

An amount of money P is invested in an account where interest is compounded at the end of a period. The future worth F yielded at an interest rate i after n periods may be determined from the following formula:F=P(1+i)^n. Write a program that will calculate the future worth of an investment for each year from 1 to n. The input to the function should include the initial investment P, the interest rate i(as a decimal), and the number of years n for which the future worth is to be calculated. Run for  P = 100,000, i = .06 and n = 5 years.

P = 100000;

i = .06;

n = 1:5; % years

F=P*(1+i).^n;

plot(F);

This produces the result of

F =

   1.0e+05 *

    1.0600    1.1236    1.1910    1.2625    1.3382

And so the answer, after 5 years, is 133820

2.7

The "divide and average" method, an old-time method for approximating the square root of any positive number a can be formulated as x=(x+a/x)/2.

For this problem, we wrote a function

function [ sqrroot ] = fp27( a )

tol = 1.0E-5;

if (a>0)

    

     x = a/2; 

    while(1==1)

      y = (x+a/x)/2;

      e = abs((y-x)/y);

      x=y;

      if (e<tol) break;

    end

    end

    

    sqrroot = x;

else

    sqrroot = 0;

end

end

which can be implemented for verification as follows

EDU>> fp27(4)

ans =

     2

EDU>> fp27(9)

ans =

    3.0000

EDU>> fp27(16)

ans =

    4.0000

2.4

The cosine function can be evaluated by the following infinite series : cos(x) = 1 - (x^2)/2!+(x^4)/4!-(x^6)/6!+ (x^8)/8! ... Write an algorithm to implement this formula.

N=50;

x=pi;  % Assumes x = pi for ease of verification

terms = [0:2:N];

facts = factorial(terms);

Power_Of_X = x.^terms;

signs = (-1).^(terms/2);

quotes = Power_Of_X ./facts;

series = signs .* quotes;

M_cos_N = cumsum(series);

plot (M_cos_N);

The problem then asks for the percent relative error or the algorithm. This can be done by updating the above code to include on the end:

%Percent Error

true = -1.0;

P_error = (true - M_cos_N)/true * 100;

figure(2);

plot(terms,(abs(P_error)));

This produces the graph:

Showing that the error inherent in the algorithm becomes largely negligible around 30 terms in

1.8

1.11

According to the problem description, 1L is ingested as food, 0.3L are produced metabolically, 0.05L inhaled, 0.4L exhaled, 0.2L lost to sweat, 1.4L lost to urine, 0.2L lost to feces, and 0.35L lost through the skin. Simplified, this is represented as -1 - 0.05 + 0.4 + 0.2 + 1.4 + 0.2 + 0.35 liters lost each day, which comes to 1.5L, meaning one must drink 1.5L of water each day to maintain steady-state condition