Embedded Files
Problem 1
Problem 1
Problem Description
Problem Description
Figure 1: A falling sphere [2]
We have gone over a number of engineering problem-solving scenarios involving an object flowing through a given fluid. This fluid can be any gas or liquid medium. In the case of the falling parachutist, for example, the person represents the object while the air represents the liquid. An object, as you can imagine, performs a level of disturbance to the fluid as it falls through it. The measure of this disturbance can be classified into conditions. Two that we will explore are laminar and turbulent flow. In order to mathematically determine the disturbance of the object falling through the liquid, we use a value called Reynolds number, Re. This calculation as outlined below considers the density of the fluid, ρ, the dynamic viscosity of the fluid, μ, the diameter of the object, d, and the object's settling velocity, v.
When Re < 0.1, the falling object is said to be under laminar conditions. Alternatively, when Re >= 0.1, it is under turbulent conditions. To find the settling velocity of the object, we must first consider the Reynolds number. Based on the conditions, we can use one of two equations. The first below is derived from Stokes Law and is used for laminar conditions. In addition to the values outlined above, this equation uses the object's density, ρs and the gravitational constant, g.
When the conditions are turbulent, we use the following equation which uses the drag coefficient, CD.
Just as the conditions are determined based on Reynolds number, the drag coefficient uses equation 1. The equation for CD is as follows.
In this problem solving lab, we use the information outlined above along with the following conditions to find the velocity of a falling particle.
Table 1: Initial Conditions
Numerical methods are needed to solve this problem because (as will be shown in the alternative approach), the complexity of the equations does not allow the velocity to simply be found. An estimate can be found using the Stokes equation, but more accurate numerical methods are necessary to solve for the velocity. Furthermore, if Stokes equation is used to estimate velocity and this estimate is, in turn, used to find the Reynolds number, it will be found that flow is turbulent flow. However, the Stokes equation assumes that flow is laminar. This is further indication that a numerical method is needed to solve for the velocity.
Numerical Approach
Numerical Approach
Fixed-Point Iteration
Fixed-Point Iteration
Modified Secant Method
Modified Secant Method
The first approach used to solve this problem was fixed-point iteration. In this method, the goal is to identify two equations, denoted f(v) and g(v), that can be used to estimate the root. A new value of the root is estimated based on the value of g at the old estimate. In other words, a new estimate is found by evaluating g at the old estimate, or:
Like all of the employed methods discussed on this page, the process is repeated until the root is designated to be within a certain error criterion.
Sample Calculation
Sample Calculation
The following is a sample calculation for the fixed-point method. The initial guess was v = 0.5. With this method,
Using the known values for the constants and the initial guess of 0.5 m/s:
The error at this new point can then be evaluated to see if another iteration is needed.
The goal of the modified secant method is to use a fractional perpetuation of the independent variable (v) to estimate the derivative of the function. The fractional perpetuation is denoted by the symbol delta. This gives the equation:
This f'(v) function can be applied to an equation for the Newton-Raphson method (a different method not employed on this page), and with some algebraic manipulation, it can be said that:
This estimate for the next x value is how the modified secant method is emplyed until the stopping criterion is reached.
Sample Calculation
Sample Calculation
The following is a sample calculation for the modified secant method. The d or step size, was defined to be 0.01. The initial guess was chosen to be 0.5 m/s.
The error at this new point can then be evaluated to see if another iteration is needed.
Click for the numerical approach pseudo-code using the fixed point iteration method
Click for the numerical approach pseudo-code using the fixed point iteration method
%% Declare variables
delta = 0.01; % Delta value
vr = 0.04; % Initial guess based on plot
es = 0.0001; % Stopping error
iter = 0; % Iteration value
ea = 100; % Approximate error
maxit = 100; % Maximum number of iterations
d = 0.00035; % Distance (m)
rho = 1000; % Density of liquid (kg /m^3)
rhoS = 7874; % Density of object (kg /m^3)
mu = 0.0014; % Dynamic viscosity (kg / m*s)
g = 9.81; % Gravitational constant (m/s^2)
h = 0.001; % Step size
total = 0.2; % Total iterations for vector v
FUNCTION fixedPointMeth(vr, es, iter, ea, maxit, d, rho, rhoS, mu, g)
DO
IF (iter < maxit) && (ea >= es)
THEN vrold = vr;
vr = sqrt(((4.*g.*(rhoS-rho).*d)./(3.*rho.*(((24.*mu)./(rho.*vrold.*d))+(3./sqrt(rho.*vrold.*d./mu)))+(0.34))));
IF(vr ~= 0)
THEN ea = (abs(vr - vrold)/vr);
ELSE break
END IF
iter = iter + 1;
END IF
fvr = sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*vr.*d)./mu)) + (3./sqrt(((rho.*vr.*d)./mu))) + 0.34).*rho)) - vr;
root = vr;
END DO
END FUNCTION
Click for the numerical approach pseudo-code using the modified secant method
Click for the numerical approach pseudo-code using the modified secant method
%% Declare variables
delta = 0.01; % Delta value
vr = 0.04; % Initial guess based on plot
es = 0.0001; % Stopping error
iter = 0; % Iteration value
ea = 100; % Approximate error
maxit = 100; % Maximum number of iterations
d = 0.00035; % Distance (m)
rho = 1000; % Density of liquid (kg /m^3)
rhoS = 7874; % Density of object (kg /m^3)
mu = 0.0014; % Dynamic viscosity (kg / m*s)
g = 9.81; % Gravitational constant (m/s^2)
h = 0.001; % Step size
total = 0.2; % Total iterations for vector v
FUNCTION secantMeth(delta, vr, es, iter, ea, maxit, d, rho, rhoS, mu, g)
DO
IF (iter < maxit && ea >= es)
THEN vrold = vr;
vr = vrold - ((delta).*(vrold).*(sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*vrold.*d)./mu)) + (3./sqrt(((rho.*vrold.*d)./mu))) + 0.34).*rho)) - vrold))./((sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*(vrold + (delta.*vrold)).*d)./mu)) + (3./sqrt(((rho.*(vrold + (delta.*vrold)).*d)./mu))) + 0.34).*rho)) - (vrold + (delta.*vrold)))-(sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*(vrold).*d)./mu)) + (3./sqrt(((rho.*(vrold).*d)./mu))) + 0.34).*rho)) - (vrold)));
IF(vr ~= 0)
THEN ea = abs((vr - vrold) /vr);
ELSE
ea = 0;
END IF
iter = iter + 1;
END IF
fvr = sqrt((4*g*(rhoS - rho)*d)/(3*((24/((rho*vr*d)/mu)) + (3/sqrt(((rho*vr*d)/mu))) + 0.34)*rho)) - vr;
root = vr;
END DO
END FUNCTION
Alternative Approach
Alternative Approach
Our goal in this problem solving lab is to find the roots of the velocity equation to determine the values of the settling velocity of a falling particle under the stated conditions. To do so, we can begin by taking the following three conditions and combining them to find a value in terms of v.
Combining these equations into one equation:
Since
Plugging CD into the equation for v then simplifying yields
We can make an initial value of the root to get us started. To do so, we can use equation 2 and the assumed values to make the following calculation.
From this value, we can conclude that the root falls somewhere between 0 and 1. With this information, we employ the following pseudo-code to graph the function.
Click for the alternative approach pseudo-code
Click for the alternative approach pseudo-code
d = 0.00035
rho = 1000
rhoS = 7874
mu = 0.0014
g = 9.81
h = 0.001
total = 10
iter = 0
Re = zeros(total, 1)
CD = zeros(total, 1)
y = zeros(total, 1)
v = zeros(total, 1)
DO
IF (iter ~= 1) THEN v(iter, 1) = v(iter - 1, 1) + h
END
Re(iter, 1) = (rho*v(iter, 1)*d)/mu
CD(iter, 1) = (24/Re(iter, 1)) + (3/sqrt(Re(iter, 1))) + 0.34
y(iter, 1) = sqrt((4*g*(rhoS - rho)*d)/(3*CD(iter, 1)*rho)) - v(iter, 1)
iter = iter + 1
IF (iter == 1000) THEN EXIT
END IF
END DO
plot(v, y)
set(gca, 'FontSize', 14)
xlabel('Velocity (m/s)')
ylabel('f(v)')
grid on
Results
Results
Alternative Solution Results
Alternative Solution Results
From the alternative solution steps and pseudo-code, we obtain the following plots.
Figure 2: The function of velocity, v
Figure 3: The function of velocity, v, Zoomed-in
After zooming in, we can see that the estimate is v ≈ 0.14475037608 m/s.
Numerical Solution Results
Numerical Solution Results
The figure below shows the function of the velocity with marks that represent the location of estimates using open methods.
Figure 4: The function graphed with the calculated roots using the fixed-point and modifed secant methods.
Figure 5: The function graphed with the calculated roots using the fixed-point and modifed secant methods - zoomed in to see where the function crosses the x-axis
The fixed-point method estimates the root to be 0.1447455 m/s while the modified secant method estimates the root to be 0.14475065 m/s.
Error Analysis
Error Analysis
As stated above, the alternative method estimates the root to be about 0.1447503760 m/s. If we assume this value to be the true value, the following true error calculations can be made.
Fixed-Point Method
Fixed-Point Method
Modified Secant Method
Modified Secant Method
As seen above, the modified secant method is more accurate, as its true error is lower than that of the fixed-point iteration method. As we will see now, the number of iteratons and relative errors of the methods also suggest that the modified secant method is the best estimate out of these two open methods.
Figure 6: The relative error as the iteration increases for fixed-point method
Figure 7: The relative error as the iteration increases for modified secant method
The above graphs provide information about the convergence of the two open numerical methods. The convergence of the fixed point method was affected by the intial guess of the root. As can be seen from the plot on the left, relative error decreased with each of the 8 iterations. Even though it required more iterations, the method did converge to a point.
The convergence of the modified secant method was faster, as only 4 iterations were required. Figure 7 shows that this method oscillated around the root before it converged (evidenced by the rise in relative error at the second iteration).
As seen in the pseudocode, each estimate estimates the root until the relative error is calculated to be less than some chosen error criterion. In this case, the stopping error chosen, es, was initiated to a value of 0.0001. In other words, the method would continue to run iterations until this error condition was met. The fixed-point method ran through 8 iterations before the relative error was found to be 8.0978e-05 while the modified secant method ran through 4 iterations before stopping at a relative error of 2.4508*10^(-5). From these values, it can be concluded that the modified secant method has less approximation error while also running for less time. On the other hand, the fixed-point method runs for more iterations and still has more approximate error when it stops running. It is important to note the way in which the relative error was calculated. Since we have no real value for the root, otherwise we wouldn't have to use these estimations, we take the relative error as a relationship between the current estimation value of velocity and the previous one as summarized in the following equation.
An example calculation of the relative error is as follows.
The approximate error is above the error threshold of 0.1%. Therefore, a new estimate for the root would be found using one of the above equations, depending on the open method being performed. This would continue until the approximate error was less than 0.1% (or the designated stopping criterion).
Discussion
Discussion
Eight iterations of the fixed-point method and four iterations of the modified secant method were required to meet the error criterion of 0.0001. As explained above, the modified secant method also had less relative error associated with the final root value. Therefore, it can be concluded that the modified secant better was the better approximation for the root of the derived equation.
Instead of directly solving for velocity using Stokes law, a numerical method was used to minimize error and get as close to the true value as possible. The Stokes equation estimate the root to be 0.3278 m/s, but the true value of the root was 0.144750 m/s. The emplyed numerical methods were able to find the true root within 0.001953%. It is clear that the numerical method (specifically the modified secant method) was the better approach because of its accuracy and relatively small number of iterations.
The rate of convergence of these methods is dependent of the initial guesses chosen because the initial guess is used in the equations for both methods. In order to get the next estimate (that is theoretically closer to the true value), the initial guess must be used. To minimize the number of iterations needed, the initial estimate should be as close to the true value as possible. This is why the estimate based on Stokes law is useful. The stopping criterion also affects the rate of convergence. If a larger, less precise stopping criterion is chosen, then fewer iterations are needed because it allows the program a wider range of acceptable roots. To get a more accurate and precise result, the initial estimate and the stopping criterion should be carefully chosen.
Click here for the full MATLAB code used to solve Problem 1
Click here for the full MATLAB code used to solve Problem 1
%% Declare variables
delta = 0.01; % Delta value
vr = 0.04; % Initial guess based on plot
es = 0.0001; % Stopping error
iter = 0; % Iteration value
ea = 100; % Approximate error
maxit = 100; % Maximum number of iterations
d = 0.00035; % Distance (m)
rho = 1000; % Density of liquid (kg /m^3)
rhoS = 7874; % Density of object (kg /m^3)
mu = 0.0014; % Dynamic viscosity (kg / m*s)
g = 9.81; % Gravitational constant (m/s^2)
h = 0.001; % Step size
total = 0.2; % Total iterations for vector v
%% Declare Equations
v = (0:h:total);
f = sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*v.*d)./mu)) + (3./sqrt(((rho.*v.*d)./mu))) + 0.34).*rho)) - v;
%% Plot Equations
figure(1)
plot(v, f, '-b', 'LineWidth', 1);
title('Velocity function');
ylabel('f(v)');
xlabel('Velocity (m/s)');
grid on
hold on
plot(0.14475037608, 0, 'r*');
title('Velocity Function with Assumption of the Root');
legend('Function', 'Assumption of the Root');
set(gca, 'FontSize', 14);
figure(2)
plot(v, f, '-b', 'LineWidth', 1);
title('Velocity function');
ylabel('f(v)');
xlabel('Velocity (m/s)');
set(gca, 'FontSize', 14);
grid on
hold on
%% Function Calls
fixedPointMeth(vr, es, iter, ea, maxit, d, rho, rhoS, mu, g);
secantMeth(delta, vr, es, iter, ea, maxit, d, rho, rhoS, mu, g);
%% Implement fixed-point iteration
function fixedPointMeth(vr, es, iter, ea, maxit, d, rho, rhoS, mu, g)
eaVector = zeros(8, 0);
iterVector = zeros(8, 0);
while (iter < maxit) && (ea >= es)
vrold = vr;
vr = sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*vrold.*d)./mu)) + (3./sqrt(((rho.*vrold.*d)./mu))) + 0.34).*rho));
eaVector(iter + 1, 1) = ea;
iterVector(iter + 1, 1) = iter;
iter = iter + 1;
if (vr ~= 0)
ea = (abs(vr - vrold)/vr);
else
break
end
end
figure(4)
plot(iterVector, eaVector, 'b', 'LineWidth', 1);
hold on
plot(iter - 1, ea, 'rx', 'LineWidth', 1);
title('Approximate Relative Error as Iterations Increase for Fixed-Point Method')
xlabel('Number of Iterations');
ylabel('Approximate Relative Error');
legend('Approximate Relative Error', 'Final Approximate Relative Error');
grid on
xlim([0 9]);
set(gca, 'FontSize', 14);
figure(2)
fvr = sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*vr.*d)./mu)) + (3./sqrt(((rho.*vr.*d)./mu))) + 0.34).*rho)) - vr;
root = vr;
format long
plot(vr, fvr, 'k*');
hold on
disp(root);
disp(ea);
end
%% Implement Modified Secant Method
function secantMeth(delta, vr, es, iter, ea, maxit, d, rho, rhoS, mu, g)
eaVector = zeros(4, 0);
iterVector = zeros(4, 0);
while (iter < maxit && ea >= es)
vrold = vr;
vr = vrold - ((delta).*(vrold).*(sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*vrold.*d)./mu)) + (3./sqrt(((rho.*vrold.*d)./mu))) + 0.34).*rho)) - vrold))./((sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*(vrold + (delta.*vrold)).*d)./mu)) + (3./sqrt(((rho.*(vrold + (delta.*vrold)).*d)./mu))) + 0.34).*rho)) - (vrold + (delta.*vrold)))-(sqrt((4.*g.*(rhoS - rho).*d)./(3.*((24./((rho.*(vrold).*d)./mu)) + (3./sqrt(((rho.*(vrold).*d)./mu))) + 0.34).*rho)) - (vrold)));
eaVector(iter + 1, 1) = ea;
iterVector(iter + 1, 1) = iter;
iter = iter + 1;
if (vr ~= 0)
ea = abs((vr - vrold) /vr);
else
ea = 0;
end
end
figure(3)
plot(iterVector, eaVector, 'b', 'LineWidth', 1);
hold on
plot(iter - 1, ea, 'rx', 'LineWidth', 1);
title('Approximate Relative Error as Iterations Increase for Modified Secant Method')
xlabel('Number of Iterations');
ylabel('Approximate Relative Error');
legend('Approximate Relative Error', 'Final Approximate Relative Error');
xlim([0 5]);
grid on
set(gca, 'FontSize', 14);
figure(2)
fvr = sqrt((4*g*(rhoS - rho)*d)/(3*((24/((rho*vr*d)/mu)) + (3/sqrt(((rho*vr*d)/mu))) + 0.34)*rho)) - vr;
root = vr;
disp(root);
format long
plot(vr, fvr, 'r*');
legend('Function', 'Fixed-Point Estimate', 'Modified Secant Estimate');
title('Velocity Function with Root Estimates using Open Methods');
xlabel('Velocity (m/s)');
ylabel('Function of Velocity');
set(gca, 'FontSize', 14);
disp(root);
disp(ea);
end
Problem 2
Problem 2
Problem Description
Problem Description
In engineering 311, we worked with RLC circuits which uses three components: a resistor, an inductor, and a capacitor. the figure below shows this RLC circuit. [1]
Figure 8: An RLC circuit
Impedance is used to describe the resistance of the components with resistances that depend on angular frequency. That being said, the following equation derived from Kirchhoff's laws is used to calculate the impedance, Z, based on resistance, R, capacitance, C, inductance, L, and the angular frequency, ω.
The following table shows the initial conditions of an RLC circuit.
Table 2: Initial Conditions
The task of this problem is to find the ω value in which Z = 75 Ω. To do so, we implement two types of solutions. The first, the numerical solution, investigates the bisection and false-position bracketing methods. An alternative solution looks at the graph of the function to estimate the ω value when Z = 75 Ω.
Numerical Approach
Numerical Approach
As stated in the problem description, we are given the following values:
Z = 75 Ω
xl = 1 radians/second
xu = 1000 radians/second
R = 225 Ω
C = 0.6 * 10^{-6} F
L = 0.5 H
We are given the following equation:
Now, we can display this equation to be in the form f(x) = 0 as follows
Now that our equation is in the form f(ω) = 0, we can implement two different bracketing methods to find the root of this equation.
Bisection Method
Bisection Method
The first method we explore is the bisection method. For this method, two values, xl and xu are chosen to bracket the root. As outlined in the problem statement, the initial bounds chosen are 1 radian/second and 1000 radians/second, respectively. We begin this method by finding the value of f(x) when ω = xl and ω = xu. Then, we find the estimate of the root, xr by using the following equation.
In the next step, the value of f(xr) is determined. If this value is greater than 0, then we move the upper bound to the estimated root and the lower bound stays the same. Alternatively, if the value is less than 0, we move the lower bound to the estimated root and the upper bound stays the same. Lastly, if the value is 0, this means we found the root and we need no further calculations. When the function value at xr is not equal to 0, we cycle back through this procedure, changing the bound values each time. We keep a counter for the number of iterations called iter. This counter is updated each time the code cycles through, which keeps track of the amount of times the calculations were done. Also in each iteration, we calculate an error value. Since we don't know the actual value of the root, our error equation uses the current and previous values of xr. The following equation is used to find the relative error.
If the ea falls below a predetermined error criterion, we can assume that the root is found to some level of accuracy, depending on the criterion chosen.
Sample Calculation
Sample Calculation
Because -0.02 is a negative number, it is known that there is a root within the bounds of 1 and 1000. Next, a new estimate for x is calculated.
Now, the following calculation is done and because it yields a negative value, the upper bound is replaced with the new x estimate.
The new bounds are 1 and 500.5 rad/s.
False-Position Method
False-Position Method
The false-position method follows the same procedure as the bisection method with the exception of the calculation for xr as the iteration increases. For the false-position method, the following equation is used.
Other than the change to this equation, the false-position method solution steps follow those of the bisection method.
Sample Calculation
Sample Calculation
The following is a sample calculation for the false position method, starting with the initial bounds of 1 and 1000 rad/s. The function is evaluated at 1 and 1000 and yields the following values.
The new x estimate was found using the following formula.
The function evaluated at this new x value yields -0.0087 rad/s. Because this is a negative value, the upper limit is replaced with this new x. The bounds become 1 and 995.6574 rad/s.
Click here for the Numerical Solution Pseudo-code using Bisection Method
Click here for the Numerical Solution Pseudo-code using Bisection Method
R = 225
C = 0.6*(10^(-6))
L = 0.5
xl = 1
xu = 1000
es = 0.001
imax = 1000
xr = 3000
iter = 0
ea = 100 % ea is relative approximate error
Z = 75
bisectionMethod(xl, xu, es, imax, xr, ea, iter, R, C, L, Z);
FUNCTION bisectionMethod(xl, xu, es, imax, xr, ea, iter, R, C, L, Z)
iter = 0
fl = f(xl)
DO
xrold = xr
xr = (xl+xu)/2
fr = f(xr)
iter = iter + 1
IF (xr ~= 0) THEN ea = abs((xr - xrold) / xr) * 100
END IF
IF (fr > 0) THEN xu = xr
ELSEIF (fr < 0) THEN xl = xr
ELSE ea = 0
END IF
IF (ea < es || iter >= imax) EXIT
END IF
END DO
END FUNCTION
Click here for the Numerical Solution Pseudo-code using False-Position Method
Click here for the Numerical Solution Pseudo-code using False-Position Method
R = 225
C = 0.6*(10^(-6))
L = 0.5
xl = 1
xu = 1000
es = 0.001
imax = 1000
xr = 3000
iter = 0
ea = 100 % ea is relative approximate error
Z = 75
fPositionMethod(xl, xu, es, imax, xr, iter, ea, R, C, L, Z);
FUNCTION fPositionMethod(xl, xu, es, imax, xr, iter, ea, R, C, L, Z)
iter = 0
fl = f(xl)
fu = f(xu)
DO
xrold = xr
xr = xu - ((fu * (xl - xu))/(fl - fu))
fr = f(xr)
iter = iter + 1
IF (xr ~= 0) THEN ea = abs((xr - xrold) / xr) * 100
END IF
IF (fr > 0) THEN xu = xr
ELSEIF (fr < 0) THEN xl = xr
ELSE ea = 0
END IF
IF ea < es || iter >= imax EXIT
END IF
END DO
END FUNCTION
Alternative Approach
Alternative Approach
As an alternative solution, we can graph the function we solved for above. We can use the MATLAB pseudo-code below. In MATLAB, we can find the value of the function when it crosses the x-axis using either the data cursor or the fzero function (see Alternative Solution Results for graph). The function was found to cross the x-axis at 157.9088807 rad/s.
Click here for the Alternative Solution Pseudo-code
Click here for the Alternative Solution Pseudo-code
FOR i = 1:1:100000
IF (i == 1) THEN x(i, 1) = 0
ELSE x(i, 1) = x(i - 1, 1) + 0.1
END IF
f(i, 1) = (1 / Z) - (sqrt((1 / (R^2)) + ((x(i, 1)*C) - (1 / (x(i, 1)*L)))^2));
END FOR
PLOT (x, f, 'LineWidth', 1.5)
label x
label y
legend('Function', 'Bisection Estimate', 'False-Postition Estimate')
Results
Results
Numerical Solution Results
Numerical Solution Results
The figure below shows the function and the estimated roots from the bisection method and the false-position method.
Figure 9: Impedance vs. ω Plot
Figure 10: Impedance vs. ω Plot - Zoomed in
The figure to the right is the same is the one on the left, but the x-axis is made to be much smaller in its range so that the difference between the approximations can be observed. The bisection method approximates the velocity to be 157.9017 rad/s while the false-position method approximates the velocity to be 158.5284 rad/s.
Alternative Solution Results
Alternative Solution Results
As an alternative solution, we can graph the function and observe where it crosses the x-axis.
Figure 11: The function of velocity when it crosses the x-axis
Figure 12: The function of velocity when it crosses the x-axis, zoomed-in
As seen above, the assumption we make for the root of the function is ω ≈ 157.9088807 rad/s.
Error Analysis
Error Analysis
The result of the bisection method is more accurate than the false-position method. This conclusion can be reached from the alternative method data because the function is seen to cross at 157.9088807 radians / second, as described in the results section of the alternative solution. If we assume this to be the true value, then we can see how far the results of the bracketing methods deviate from the true value. We can use the following equation.
Bisection Method
Bisection Method
False-Position Method
False-Position Method
From the true errors, it is evident that the bisection method is a better approximation of the root than the false-position method when compared to the alternative solution estimation.
As seen in the pseudo-codes outlined above, while loops are used to guide the code. In order to ensure that precision requirements are met, one of the conditions that causes the code to continue to run is that the relative error must be greater than the predetermined error criterion. The relative error uses equation 6 above. More specifically, this equation depicts the approximate percent relative error which we can call ea. Only when ea falls below the error criterion can the code break out of this while loop. The figures below show the relative approximate error values for every iteration value for each method.
Figure 13: The False-Position Method: relative error vs. number of iterations
Figure 14: The Bisection Method: relative error vs. number of iterations
The stopping error criterion for each of these plots was 0.001, or 0.1%. The false-position method ran for 427 iterations before stopping at a relative approximate error value of 1.8696e-05 which is marked as a red 'x' on figure 13 above. On the other hand, the bisection method runs through 12 iterations before stopping at a relative approximate error value of 7.7208e-04, also marked as a red 'x' on figure 14.
The patterns of both of these error plots suggest that the methods converge to one value. Because these plots represent the relative approximate error, this value isn't necessarily the real root of the function. As outlined in the results section, the false-position method converges to 158.5283 rad/s while the bisection method converges to 157.9017 rad/s. The observation that the bisection method runs through far less iterations than the false-position method suggests that the bisection method is more time efficient.
As seen above by the number of iterations required to reach the error criterion, the bisection method had a faster rate of convergence than the false position method. This is due to the shape of the original function. The function is nonlinear, and the false-position method, because it relies on the slope of the function at a point, was stuck in an area of the graph that prevented it from quickly finding the root. This is why Figure 13 shows a long (nearly) horizontal line; the x values were not changing much in relation to each other, but they were not adequately close to the root. The shape of the graph and its effects on these methods will be further discussed in the next section.
Discussion
Discussion
12 iterations of the bisection method were required to determine the answer precisely according to the error criterion. Alternatively, the false-position method requires 426 iterations. The final error values are also outlined in this table. The number of iterations for the methods are different because they use different equations to find xr in each iteration. While the bisection method uses the average of the bounds, the false-position method approximates the slope of the line and uses this slope to find a value for xr. There are more iterations required for the false position method to approximate the root because the xr value does not change as dramatically in every iteration as compared to the bisection method. Both of these bracketing methods converge since the approximate value gets increasingly closer to the true value. The rate of this convergence is different between the two bracketing methods used in this example. The convergence rate of the bisection method is much faster than that of the false-position method as seen by the difference in iterations needed to acquire the desired solution to the necessary precision.
Further, the bisection method approximation is much more accurate, as its true error is smaller than that of the false-position method. In conclusion, the bisection method is the best bracketing method of those observed in this problem.
In comparison to the open methods used in problem 1, the bisection method did not have the most efficient rate of convergence. The modified secant method converged fastest (requiring only 4 iterations). The false position method remained the most inefficient because the shape of the function hindered its ability to quickly find the root. Open methods proved to be more efficient in this context because the function converged and only one initial guess was needed. When deciding which method to use, the shape of the function and the possibilty that open methods may diverge should be considered.
Click here for the full MATLAB code used to solve Problem 2
Click here for the full MATLAB code used to solve Problem 2
%% Declare variables
R = 225; % Resistance (Ohms)
C = 0.6*(10^(-6)); % Capacitance (F)
L = 0.5; % Inductance (H)
xl = 1; % Initial lower bound guess
xu = 1000; % Initial upper bound guess
es = 0.001; % Stopping criterion
imax = 1000; % Maximum iteration value
xr = 3000; % Middle value
iter = 0; % Iteration value
ea = 100; % Starting error
Z = 75; % Impedance value
f = zeros(1000, 1); % Initialize vector for function values
x = zeros(1000, 1); % Initialize vector for x values
%% Function calls
bisectionMethod(xl, xu, es, imax, xr, ea, iter, R, C, L, f, x, Z);
fPositionMethod(xl, xu, es, imax, xr, iter, ea, R, C, L, f, x, Z);
%% Implementing Bisection Method
function bisectionMethod(xl, xu, es, imax, xr, ea, iter, R, C, L, f, x, Z)
fl = (1 / Z) - (sqrt((1 / (R^2)) + ((xl*C) - (1 / (xl*L)))^2));
eaVector = zeros(13, 1);
iterVector = zeros(13, 1);
while (ea >= es && iter < imax)
xrold = xr;
xr = (xl+xu)/2;
fr = (1 / Z) - (sqrt((1 / (R^2)) + ((xr*C) - (1 / (xr*L)))^2));
eaVector(iter + 1, 1) = ea;
iterVector(iter + 1, 1) = iter;
iter = iter + 1;
if (xr ~= 0)
ea = abs((xr - xrold) / xr);
end
if (fr > 0)
xu = xr;
elseif (fr < 0)
xl = xr;
else
ea = 0;
end
end
disp(ea)
%% Plot error
figure(2)
plot(iterVector, eaVector - 1, 'b', 'LineWidth', 1);
title('Bisection Method Approximate Relative Error per Iteration Value');
xlabel('Number of iterations');
ylabel('Relative Approximate Error');
set(gca, 'FontSize', 14);
hold on
plot(max(iterVector), min(eaVector - 1), 'rx');
ylim([-1 100]);
xlim([-.01 13]);
legend('Approximate Relative Error for Corresponding Iterations', 'Bisection Method Estimation');
%% Plot bisection method
for i = 1:1:100000
if i == 1
x(i, 1) = 0;
else
x(i, 1) = x(i - 1, 1) + 0.1;
end
f(i, 1) = (1 / Z) - (sqrt((1 / (R^2)) + ((x(i, 1)*C) - (1 / (x(i, 1)*L)))^2));
end
figure(1)
plot(x, f, '-b', 'LineWidth', 1);
hold on
plot(xr, fr, 'r*', 'LineWidth', 1);
grid on
title('Impedance vs. Angular Velocity with Root Estimations');
set(gca, 'FontSize', 14);
end
%% Implementing False Position Method
function fPositionMethod(xl, xu, es, imax, xr, iter, ea, R, C, L, f, x, Z)
eaVector = zeros(427, 1);
iterVector = zeros(427, 1);
fl = (1 / Z) - (sqrt((1 / (R^2)) + ((xl*C) - (1 / (xl*L)))^2));
fu = (1 / Z) - (sqrt((1 / (R^2)) + ((xu*C) - (1 / (xu*L)))^2));
while (ea >= es && iter < imax)
xrold = xr;
xr = xu - ((fu * (xl - xu))/(fl - fu));
fr = (1 / Z) - (sqrt((1 / (R^2)) + ((xr*C) - (1 / (xr*L)))^2));
eaVector(iter + 1, 1) = ea;
iterVector(iter + 1, 1) = iter;
iter = iter + 1;
if (xr ~= 0)
ea = abs((xr - xrold) / xr);
end
if (fr > 0)
xu = xr;
elseif (fr < 0)
xl = xr;
else
ea = 0;
end
end
disp(ea)
%% Plot false-position method on same plot
for i = 1:1:100000
if i == 1
x(i, 1) = 0;
else
x(i, 1) = x(i - 1, 1) + 0.1;
end
f(i, 1) = (1 / Z) - (sqrt((1 / (R^2)) + ((x(i, 1)*C) - (1 / (x(i, 1)*L)))^2));
end
figure(1)
hold on
plot(xr, fr, 'k*', 'LineWidth', 1);
legend('Function', 'Bisection Estimate', 'False-Postition Estimate');
xlim([-20 1120]);
ylim([-20 5]);
grid on
xlabel('Angular Frequency (Rad/s)');
ylabel('Impedance (\Omega)');
title('Impedance vs. Angular Velocity');
set(gca, 'FontSize', 14);
figure(4)
plot(x, f, '-b', 'LineWidth', 1);
hold on
plot(157.9088807, 0, 'r*');
xlim([-20 1120]);
ylim([-20 5]);
xlabel('Angular Frequency (Rad/s)');
ylabel('Impedance (\Omega)');
title('Impedance vs. Angular Velocity Approximating when Impedance is 0 \Omega');
legend('Function', 'Assumption of Real Root');
set(gca, 'FontSize', 14);
%% Plot error
figure(3)
plot(iterVector, eaVector, 'b', 'LineWidth', 1);
title('False-position Approximate Relative Error per Iteration Value');
xlabel('Number of iterations');
ylabel('Relative Approximate Error');
set(gca, 'FontSize', 14);
hold on
plot(max(iterVector), ea, 'rx');
ylim([ea 105]);
xlim([-10 450]);
grid on
legend('Approximate Relative Error for Corresponding Iterations', 'False-Position Method Estimation');
end
Reflection and Statement of Contribution
Reflection and Statement of Contribution
Katherine Reiss and Libby Welborn collaborated to make this webpage. Parts from each person's intial attempt at the problem were added to the final webpage. Feedback from each other and from instructors lead to more detailed graphs (in the numerical and alternative solution sections) and improved MATLAB code. For example, Katherine was able to correctly code the modified secant method, while Libby was able to code the fixed-point method. Ultimately, these sections were combined into one functional code. Team members received and incorporated feedback and communicated to each other about the project's progress so that both members contributed equally.
Katherine was responsible for setting up the numerical and alternative methods and corresponding code. She also created the graphs displayed on this page. Libby was responsible for sample calculations and the discussion. Both members worked on the overall MATLAB code, psuedo-code, the error analysis, and on the webpage as a whole.
References
References
[1] EGR312: Problem Solving Lab Assignment #2
[2] Image from https://www.comsol.com/blogs/modeling-a-sphere-falling-on-a-water-surface/Page updated
Report abuse