Numerical Methods for Engineers

Data approximation and round off error detection MATLAB

2013-02-16T15:56:00.003-08:00

Truncation error: An error that occurs when an approximate number is used to represent an exact number.

Round-off error: A error that occurs when the number of significant figures is limited, thus the entire number is not able to be represented.

True value = approximate value + error

so...

Error = true value - approximate value

so...

True error equation: E_t= true value - approximate value(This is the equation for the exact error value)

Normalized True error (Magnitude of error):

True fractional relative error = true error/true value

True percent relative error = true error/true value *100%

so...

True percent relative error equation:

ε_t = true error * 100% or ε_t = true value - approximate value * 100%
true value true value

Example:

A board for a boat is 78 cm long. The actual length need to put on the boat is 80 cm. What is the true error and percent error.

True error:

E_t= 80 cm - 78 cm
E_t= 2 cm

True percent relative error:

ε_t = 2 cm * 100% = 2.5%
         80 cm

Alternatively we often will not know the true error. In this case use the approximate error.

True percent approximate error equation:

ε_a = approximate error * 100%
      approximate value


Note :When dealing with error in a problem we should have a percent tolerance of error otherwise known as the amount of allowed error that we deem acceptable. We want our percent error to be lower than our predefined percent tolerance(ε_s). To measure this we must take the absolute value of our approximate percent error(ε_a).

The commonly used percent tolerance(n represents the number of significant figures):
ε_s =(0.5 x 10^2-n)%
Note: ε_s is also sometimes called the convergence criteria

Example:

We have to have a number with 5 significant figures to accurately approximate our value.

ε_s =(0.5 x 10^2-5)% = 0.0005%

True percent approximate error equation(iterative problems):

ε_a = current approximate value - previous approximate value * 100%
                               current approximate value

Example:
We can use the Maclaurin series(Taylor series is centered at zero) to approximate e ^xat x = .8.

e ⁸is approximately 2.22554092849......

e ^x= 1 + x + x ² + x ³+ x ⁴+ x ⁵+ x ⁶+ x ⁷+
                        2       3!     4!     5!    6!     7!       ..............
estimate 1
e ⁸= 1

ε_t =  2.22554092849 - 1 * 100% = 100% error
      2.22554092849

ε_a = 1 - 0 * 100% = 100% error
           1

estimate 2
e ⁸= 1 + .8 =1.8

ε_t =  2.22554092849 - 1.8 * 100% = 19.120786458% error
      2.22554092849

ε_a = 1.8 - 1 * 100% = 44.444444444% error
           1.8

estimate 3
e ^x= 1 + .8 + .8 ² =2.12
                          2

ε_t =  2.22554092849 - 2.12 * 100% = 4.742259607% error
      2.22554092849

ε_a = 2.12 - 1.8 * 100% = 0.15094339622% error
           2.12

estimate 4
e ^x= 1 + .8 + .8 ² + .8 ³ = 2.20533333333
2        3!

ε_t =  2.22554092849 - 2.20533333333 * 100% = 0.90798578% error
      2.22554092849

ε_a =  2.20533333333 - 2.12 * 100% = 3.834273826% error
           2.20533333333

estimate 5
e ^x= 1 + .8 + .8 ² + .8 ³ + .8 ⁴ = 2.2224
2        3! 4!

ε_t=  2.22554092849 - 2.2224 * 100%  = 0.14113101447% error
        2.22554092849

ε_a= 2.2224 - 2.20533333333  * 100% = 0.76793856506 % error
  2.2224

estimate 5
e ^x= 1 + .8 + .8 ² + .8 ³ + .8 ⁴ + .8 ⁵ = 2.22513066667
2        3! 4! 5!

ε_t=  2.22554092849 - 2.22513066667 * 100%  = 0.018434251% error
        2.22554092849

ε_a= 2.22513066667 - 2.2224  * 100% = 0.122719385 % error
2.22513066667

ε_s =(0.5 x 10^2-3)% =0.05%

xxxxxxxxxxxxxxxxxxxxxxxxxx
under construction
xxxxxxxxxxxxxxxxxxxxxxxxxx
As an Amazon Associate I earn from qualifying purchases.

Simplex Method (algorithm)

2012-05-12T14:44:00.001-07:00

Simplex Method (algorithm)

Definition: A linear programming algorithm used to find the optimal solution of a problem.
The Algorithm was mad by George Bernard Dantzig.

How to solve a problem with the Simplex Method

The following is a table of constraints for a gas company producing two types of gas

We can derive the following equations from the table

Z = Total Profit = 150 * regular + 175 * premium
Total Gas Used = 7 * regular + 11 * premium =< 77
Total Production Time = 10 * regular + 8 * premium =< 80
regular =< 9
premium =< 6

regular >= 0
premium >= 0

Approximate solution

Under construction

As an Amazon Associate I earn from qualifying purchases.

Problem 16.25 Numerical Methods

2012-05-10T08:23:00.000-07:00

Problem 16.25 Numerical Methods

The torque transmitted to an induction motor is a function of the slip between the rotation of the stator field and the rotor speed s where slip is defined as

s = n - nr/n

where n = revolutions per second of rotating stator speed and nr = rotor speed. Kirchoff's laws can be used to show that the torque (expressed in dimensionless form) and slip are related by

T = @(s) 15 * (s-s .^2) ./( (1-s) .* (4*s.^2 - 3 * s +4));

Figure P16.25 show this function. Use a numerical method to determine the slip at which the maximum torque occurs.

MATLAB Graph

Problem 16.23 Numerical Methods

2012-05-10T08:22:00.001-07:00

Problem 16.23 Numerical Methods

16.23 A total charge Q is uniformly distributed around a ring shaped
conductor with a radius a. A charge q is located at a distance x from
the center of the ring (Fig. P16.23). The force exerted on the charge by
the ring is given by

F = 1/4pie0 * qQx/(x^2 + a^2)^3/2

where e0 = 8.85 x 10^-12 C^2/(N m^2), q = Q =2 x 10^-5 C, and a = .9m.
Determine the distance x where the force is maximum.

MATLAB Code

Graph

Solution

x = 0.4025

F = 7.8820e-015

Problem 16.17 Numerical Methods

2012-05-10T08:22:00.000-07:00

Problem 16.17 Numerical Methods

The streete-Phelps model can be used to compute the dissolved oxygen concentration in a river below a point discharge of sewage (figure P16.16),

o = o_s - kd * Lo/kd +ks - ka * (e^ka *t -e^-(kd + ks)*t - Sb/ka (1 - e^-ka *t)

where 0 = dissolved oxygen concentration (mg/L), os = oxygen saturation concentration (mg/L), t = travel time (d), Lo = biochemical oxygen demand (BOD) concentration at the mixing point (mg/L)' kd = rate of decomposition of BOD (d^-1), ks = rate of settling of BOD (d^-1), ka = reaeration rate (d^-1), and Sb = sediment oxygen demand (mg/L/d).

As indicated in Fig. P16.17, Eq. (P16.17) produces an oxygen "sag" that reaches a critical minimum level Oc some travel tim tc below that point discharge. This point is called "critical because it represents the location where biota that depends on oxygen (like fish) would be the most stressed. Determine the critical travel time and concentration, given the following values:

o_s = 10; mg/L
ks = 0.06; d^-1
kd = .2; d^-1
Lo = 50; mg/L
ka = .8; d^-1
Sb = 1; mg/L/d

MATLAB Code

MATLAB Graph of approximation to minimum oxygen consentration

Solution

time =

3.5833

oxygen saturation =

-18.0722

Problem 16.4 Numerical Methods

2012-05-10T08:21:00.000-07:00

Problem 16.4 Numerical Methods

The specific growth rate of yeast that produces an antibiotic is a function of the foods concentration c,

g = 2c/ 4 + .8c +c^2 + .2c^3

As depicted in Fig. P16.4, growth goes to zero at very low concentrations due to food limitation. It also goes to zero at high concentrations due to toxicity effects. Find the value of c at which growth is a maximum.

MATLAB Code

MATLAB Plot

xl =

1.5403

ans =

0.3696

Problem 15.17 Numerical Methods

2012-05-09T08:20:00.002-07:00

Problem 15.17 Numerical Methods

The length of the longest ladder that can negotiate the corner depicted in Figure P15.17 can be determined by computing the value of theta that minimizes the following function

L (theta,alpha) = w1 /sin(theta)+ w2 /sin(pi -alpha -theta)

For the case where w1 = w2 = 2 m, use a numerical method (including software) to develop a plot of L versus a range of alpha 45 to 135 degrees.

MATLAB Code

MATLAB Plot

Problem 15.12 Numerical Methods

2012-05-09T08:20:00.001-07:00

Problem 15.12 Numerical Methods

An automobile company has two versions of the same model car for sale, a two-door coupe and the full size-four door.
(a)Graphically solve how many cars of each design should be produced to maximmize profit and what that profit is.
(b) Solve the same problem with excel

x1 = 2 door
x2 = 4 door

Tp = 13500 car1 1500 car2
production time = 15car1 + 20car2 < 8000 total production time
storage = 400 > car1 350 > car2
consumer demand = 700 x1 + 500 x2 < 240000 total cars

MATLAB Code

MATLAB Graph

Excel Total Profit equation

Excel Production Time equation

Excel Consumer Demand Equation

Excel Solver Constraints

Excel Solver Results

Problem 15.4 Numerical Methods

2012-05-09T08:20:00.000-07:00

Problem 15.4 Numerical Methods

Consider the linear programming problem:

Maximize f(x,y) = 6x + 8y

subject to

5x +2y =< 40
6x +6y =< 60
2x + 4y =< 32
x >= 0
y >= 0

Obtain the solution:
(a) Graphically
(b) Using the simplex method.
(c) Using an appropriate software package (for example, Excel).

MATLAB Code

MATLAB Graph

MATLAB Simplex solutions step by step

Constraint problem solved with excel solver

Problem 14.8 Numerical Methods

2012-05-08T11:30:00.000-07:00

Problem 14.8 Numerical Methods

Perform one iteration of the optimal gradient steepest descent method to locate the minimum of

f(x, y) = -8*x+x.^2+12*y+4*y.^2-2.*x.*y

using initial guesses x = 0 and y = 0.

MATLAB code

Surface/contour plot with approximations to lowest point

Problem 14.7 Numerical Methods

2012-05-08T11:00:00.000-07:00

Problem 14.7 Numerical Methods

Perform one iteration of the steepest ascent method to locate the maximum of
f(x,y) = 3.5x +2y + x^2 - x^4 - 2xy - y^2
using initial guesses x = 0 and y = 0. Employ bisection to find the optimal step size in the gradient search direction.

MATLAB Code

MATLAB Graph

Problem 14.6 Numerical Methods

2012-05-08T10:00:00.000-07:00

Problem 14.6 Numerical Methods

Find the minimum value of

f(x,y) = (x-3)^2 + (y-2)^2
starting at x=1 and y=1,using the steepest descent method with a stopping criterion of error = 1%. Explain your results.

MATLAB Code

Graph

Problem 14.5 Numerical Methods

2012-05-08T09:00:00.000-07:00

Problem 14.5 Numerical Methods

Find the gradient vector and Hessian matrix for each of the following functions.

(a) f(x,y) = 2xy^2 +3e^xy
(b) f = (x,y,z) x .^2 + y .^2 + 2 .*z .^2
(c) f = (x,y) log(x .^2 + 2 * x * y + 3 .*y .^2)

MATLAB Code/Soltion

Golden-Section Search

2012-05-08T02:00:00.000-07:00

Golden-Section Search
Defined as a general purpose search technique for finding a local maximum or minimum of a unimodal function

How to find the Golden Ratio

The sum of the sub lengths l₁ and l₂ equal the original length l_0.

l₀ = l_{1 +}l₂

The ratio of the lengths must be equal.

l_1/ l_{0 =}l_2/ l₁

Take the first equation and substitute it into the the denominator of the second.

l_1/ (l_{1 +}l₂)₌l_2/ l₁

Take the reciprocal of the equation.

(l_1
+l₂)_/ l₁ ₌l_1/ l₂

Now set R equal to l_2/ l_{1 .}

_R
=l_2/ l₁

Which if we apply to the entire equation will result in...

1+R = 1/R

_{Now multiply both sides of the equation by R.}

R + R² = R/R

Simplify and set the equation to zero.

R²+R -1= 0

Notice this is in the form of the quadratic equation

ax² + bx + c

Since this is in the form of a quadratic we can set it up in the form of the quadratic formula

Here is the equation set up in the quadratic formula

Solve the equation for the positive root and we arrive at the Golden Ratio

= 0.6180339887..........

As an Amazon Associate I earn from qualifying purchases.

How to use the Golden Ratio to find the local maximum

Part 1 Golden section setup

Part 2 Golden section search loop

Golden Section search Solution

Golden section search after the approximating values have been set

Golden section search after 20 iterations

Problem 13.22 Numerical Methods

2012-05-08T01:00:00.000-07:00

Problem 13.22 Numerical Methods

Use the golden-section search to determine the length of the shortest
ladder that reaches from the ground over the fence to touch the building's
wall (figure P13.22). Test it for the case where h = d = 4m

MATLAB Code

Golden Section Graph

Problem 13.21 Numerical Methods

2012-05-07T23:00:00.000-07:00

Problem 13.21 Numerical Methods

An object can be projected upward at a specified velocity. If it is subject to linear drag, its altitude as a function of time can be computed as

z = z0 + m/c(v0 _mg/c) * (1 - e^-(c/m) * t) - mg/c * t

where z = altitude (m) above the earths's surface (defined as z = 0), z0 = the initial altitude (m), m = mass(kg), c= a linear drag coefficient (kg/s), v0 = initial velocity (m/s), and t = time (s). Note that for this formulation, positive velocity is considered to be the upward direction. Given the following parameter values: g = 9.81 m/s^2, z0 = 100m, v0 = 55 m/s, m = 80kg,and c = 15 kg/s, the equation can be used to calculate the jumpers altitude. Determine the time and altitude of peak elevation (a) graphically, (b) analytically, and (c) with golden section search until the approximate error falls below 1% with initial guess of tl = 0 and tu = 10s.

Analytic Graph

Golden section Search

Solution

jumpers altitude = 192.8537 m

seconds = 3.8699

Problem 13.17 Numerical Methods

2012-05-07T20:30:00.000-07:00

Problem 13.17 Numerical Methods

The trajectory of a ball can be computed with

y = (tan(theta0))*x - g/2*v0^2*cos(theta0)^2)*x.^2 + y0;

where y = the height (m), theta0 = the initial angle(radians), v0 = the initial velocity (m/s), g = the gravitational constant = 9.81 m/s^2, and y0 = the initial height (m). Use the golden-section search to determine the maximum height given y0 = 1 m, v0 = 25 m/s and theta0 = 50 degrees. Iterate until the approximate error falls below error = 1% using initial guesses of xl = 0 and xu = 60m.

MATLAB Code

Golden Section Search

Maximum found using Parabolic Interpolation

Maximum found using Newton's Method for finding optimum

Solution

y = 19.6934

x = 31.3713

Problem 13.8 Numerical Methods

2012-05-07T08:16:00.000-07:00

Problem 13.8 Numerical Methods

Employ the following methods to find the maximum of the function from Prob 13.7

f(x) = -x^4 - 2x -8x^2 - 5x

(a) Golden section search ( xl = -2, xu = 1, error = 1%)
(b) Parabolic interpolation ( x0 = -2, xl = - 1, x2 = 1, iterations = 4). Select new points sequentially as in the secant method.
(c) Newton's method (x0 = -1, error = 1%

As an Amazon Associate I earn from qualifying purchases.

MATLAB Code Golden section search

MATLAB Code Parabolic interpolation and Newton's Method

Golden section search plot after points were set up

Golden Section search plot showing approximation to local maximum

Parabolic Interpolation

Newton's Method

Solution

Maximum = 1.4977 at x = -.4192

Problem 13.7 Numerical Methods

2012-05-06T12:27:00.000-07:00

Problem 13.7 Numerical Methods

Consider the following function:

f(x) = -x^4 - 2x -8x^2 - 5x

Use analytical and graphical methods to show the function has a maximum for some value of x in the range -2<= x <= 1

MATLAB Code

Golden section search plot after points were set up

Golden Section search plot showing approximation to local maximum

Solution

Maximum = 1.4977 at x = -.4192

Problem 12.38 Numerical Methods

2012-04-15T11:47:00.001-07:00

Problem 12.38 Numerical Methods

12.38 Linear algebraic equations can arise in the solution of differential equations. For example, the following differential equation derives from heat balance for a long, thin rod (Fig. P 12.38):

d^2T/dx^2 + h'(Ta - T) = 0

where T = temperature( degrees celsius), x = distance along the rod (m) , h' = a heat transfer coefficient between the rod and ambient air (m^-2), and T_a = the temperature of the surrounding air(C). This equation can be transformed into a set of linear algebraic equations by using a finite divided difference approximation for the second derivative (recall section 4.1.3),

d^T/dx^2 = (Ti+1 - 2Ti + Ti - 1)/ deltax^2

where Ti designates the temperature at node i. This approximation can be substituted into Eq. (P12.38) to give

-Ti - 1 + (2 + h' deltax^2)Ti - Ti+1 = h' deltax^2 Ta

This equation can be written for each of the interior nodes of the rod resulting in a tridiagonal system of equations. The first and last nodes at the rods ends are fixed by boundary conditions.

(a) Develop an analytical solution for Eq. (12.38) for a 10-m rod with Ta =
20, T(x=0) = 40, T(x = 10) = 200 and h' = 0.02.

(b) Develop a numerical solution for the same parameter values employed in (a) using a finite-difference solution with four interior nodes as shown in Fig . P 12.38( deltax = 2 m).

hp = 0.02
dx = 2;
Ta = 20

T_xe0 = 40;
T_xe10 = 200;

b = [T_xe0; ...
hp*dx^2*Ta; ...
hp*dx^2*Ta; ...
hp*dx^2*Ta; ...
hp*dx^2*Ta; ...
T_xe10]

A = [ 1 0 0 0 0 0; ...
-1 2+hp*dx^2 -1 0 0 0;
0 -1 2+hp*dx^2 -1 0 0 ;
0 0 -1 2+hp*dx^2 -1 0 ;
0 0 0 -1 2+hp*dx^2 -1 ;
0 0 0 0 0 1]

c_2 = -25.25
c_1 = 20-c_2

f_x = @(x) c_1* exp(0.141421* x)+c_2* exp(-0.141421* x)+20.
x=[0:0.01:10];
plot(x,f_x(x));

% c_1 + c_2 = 20
% c_1*exp(0.141421* 10)+ c_2* exp(-0.141421* 10)+20 = 200
% (20-c_2)*exp(0.141421* 10)+ c_2* exp(-0.141421* 10)+20 = 200

Problem 12.28

2012-04-15T11:44:00.005-07:00

%i1
%j2
%i3

%k voltage law the summation of the voltages of a closed path is zero.
% this is comprable to work like walking in a circle

%20 - 4(i2 - j2) - 2(i2 - i3) = 0

%i3: -2(i3 - i2) - 8(i3 - j2) - 5 * i3 = 0

%j2: -4(j2 - i2) - 6(j2) - 8(j2 - i3) = 0

%[i2 j2 i3]

% i1 i2 i3 i1 j2 i3
A = [ 6 -4 -2; 2 8 -15; 4 -18 8]

b = [20; 0; 0]; % from first eqaution set to 20

A\b

%20v - 6j2 - 5 i3 = 0

A =

6 -4 -2
2 8 -15
4 -18 8

ans =

5.1759
1.9095
1.7085

Problem 12.25

2012-04-15T11:44:00.003-07:00

%12.25
%Solve the circuit in Fig. P12.25 for the currents in each wire. Use gauss
%elimination with pivoting.

%Kirchhoff's current law
%At any node (junction) in an electrical circuit, the sum of currents
%flowing into that node is equal to the sum of currents flowing out of
%that node which is equal to zero.

%x1 =
%x2 = i - i2 - i3 = 0
%x3 = i3 - i4 - i5 = 0
%x4
%x5 = i2 + i7 - i8 = 0
%x6 = i4 + i6 - i7 = 0
%x7 = i5 + i9 - i6 = 0
%x8 = i8 - i9 - i10 = 0
%x9 =

%Kirchhoff's voltage law
%kvl = 110v -5i3 - 20i5 - 5i5 + 0i9 + 15i9 - 40 = 0

%kv1 = 110v -5i3 - 20i5 - 5i5 + 0i9 + 15i9 - 40v = 0 entire loop
%kvl2 = 10i7 +10i6 +0i9 +15i9 - 50i8 = 0 right loop
%kvl3 = -5i3 - 5i4 - 10i7 + 20i2 = 0 left bottom
%kvl4 = -5i5 - 5i5 - 10i6 +5i4 = 0 left top

%i1 i2 i3 i4 i5 i6 i7 i8 i9 i10
A = [0 0 -5 0 -25 0 0 0 15 -40; ...
0 0 0 0 0 10 10 -50 0 0;...
0 20 -5 -5 0 0 10 0 0 0; ...
0 0 0 5 -25 -10 0 0 0 0;];

B = [110; 0; 0; 0];

A\B

A =

6 -4 -2
2 8 -15
4 -18 8

ans =

5.1759
1.9095
1.7085

ans =

0
0
0
0
0
0
0
0
0
-2.7500

Problem 12.19 Numerical Methods

2012-04-15T11:44:00.001-07:00

Problem 12.19 Numerical Methods

Solve for the forces and reaction for the truss in Fig. P12.19. Determine the matrix inverse for the system. Does the vertical member force in the middle member seem reasonable? Why?

(Click photo to enlarge)

First we must find the sum of the forces at each point on the truss

(1)

sum of Fx = 0

F1x + F16 + F12cos(60) = 0

sum of Fy = 0

F1y + F12 * sin(60) = 0

(2)

sum of Fx = 0

F12 * cos(60) - F26 * cos(60) + F23 = 0

sum of Fy = 0

F12 * sin(60) + F26 * sin(60) = 0

(3)

sum of Fx = 0

F23 + F34 = 0

sum of Fy = 0

F36 = 0

(4)

sum of Fx = 0

F34 + F46 * cos(45) - F45 * cos(45) = 0

sum of Fy = 0

F46 * sin(45) + F45 * sin(45) = 0

(5)

sum of Fx = 0

F65 - F45 * cos(45)

sum of Fy = 0

F5y + F45 * sin(45)

(6)

sum of Fx = 0

0 = F16 + F65 + F46 * cos(45) - F26 *cos(60)

sum of Fy = 0

3500 = F26 *sin(60) + F36 + F46 * sin(45)

MATLAB Code

(Click photo to enlarge)

MATLAB solution

A = Matrix of the forces

B = Forces

C = Solution Matrix (actual forces at positions)

ans = the inverse of the A matrix

(Click photo to enlarge)

Problem 12.11 Numerical Methods

2012-04-15T11:43:00.000-07:00

Problem 12.11 Numerical Methods

A peristalic pump delivers a unit flow (Q1) of a highly viscous fluid.The network is depicted in Fig. P12.11. Every pipe section has the same length a diameter. The mass and mechanical energy balance can be simplified to obtain the flows in every pipe. Solve the following system of equations to obtain the flow in every stream.

Q3 + 2Q4 - 2Q2 = 0
Q5 + 2Q6 - 2Q4 = 0
3Q7 - 2Q6 = 0
Q1 = Q2 + Q3 == % -Q1 + Q2 + Q3 = 0
Q3 = Q4 + Q5 == % -Q3 + Q4 + Q5 = 0
Q5 = Q6 + Q7 == % -Q5 + Q6 + Q7 = 0
Q1 = Q2 + Q3

MATLAB

Solution

x = solved for the system relative to Q1

Problem 11.28 Numerical Methods

2012-04-15T10:58:00.003-07:00

Problem 11.28 Numerical Methods

A pentadiagnol system with a bandwith of five can be expressed generally as

f1 g1 h1 x1 r1
e2 f2 g2 h2 x2 r2
d3 e3 f3 g3 h3 x3 r3
... . = .
... . .
.............dn-1 en-1 fn-1 gn-1 xn-1 rn-1
..................dn en fn xn rn

Develop a program to efficiently solve such systems without pivoting in a similar fashion to the algorithm used for traditional matrices in Sec 11.11. Test it for the following case:

MATLAB Code

Solution

A = Matrix A

B = Matrix B

X = A\B

XX = linsolve(A\B)