Need help creating a function m-file for the bisection method for the following equation...

I'm really stuck on an assignment that regards the bisection method. This is due to me missing all of my lectures for the module because I have transferred course from Chemical to Mechanical Engineering. I can't get help from within the University as it is outside of term-time now and I am back home.

I am having great difficulty with the following task for which I have been given these values:

k = 1; n = 0.6833; U = 1.0714; Pgrad = -3.333333.....

I am then given the following equation defining the pressure gradient in a pipe of radius 'R' (the only unknown in the equation):

Pgrad = -2k * ( ( U * (3n + 1) ) / ( n * R ^ (1 + 1/n) ) ) ^ n

Use the bisection method to determine an approximate value of R, with an absolute error of less than 1 x 10^-6.

I'm not expecting someone to come up with the whole m-file or anything, but any help whatsoever would be greatly appreciated. I've spent hours and hours trying to catch up with where I should be on the course and even resorted to watching videos on Youtube to see if any could help.

Even if people could just try and point me in the right direction.

Thanks

(Here is the question from the actual assignment sheet)

[Rothanak]

Find the value of 3: a. Write in file M file of it define function. b. Find the root using bisection method in interval [1,2].

Answer 1

[Roger Stafford]

I recommend you carefully read the Wikipedia article at:

http://en.wikipedia.org/wiki/Bisection_method

The method requires that you first find two values a and b such that the function f(R ) whose root you are attempting to find has opposite signs for f(a) and f(b). For that purpose I advise you to make a plot of your f(R ) with varying R so you can see approximately where it crosses the zero point and thus how to set those initial starting values for a and b. For your looping until the error is less than 1e-6, you will need to use matlab's 'while' command to continue until that condition is satisfied. The above article gives a fairly clear description of how to proceed within the while-loop.

In Mathworks' File Exchange there are several contributions whose codes could possibly guide you to a solution if you get desperate. See the list at:

 http://www.mathworks.com/matlabcentral/fileexchange/index?term=bisection&utf8=✓

[Roger Stafford]

As Walter has pointed out, it is possible to solve for R explicitly, and in fact it is fairly easy to do. The manipulation goes like this:

Pgrad = -2*k*(U*(3*n+1)/n/R^(1+1/n))^n
-Pgrad/2/k = (U*(3*n+1)/n/R^(1+1/n))^n
(-Pgrad/2/k)^(1/n) = U*(3*n+1)/n/R^(1+1/n)
R^(1+1/n) = U*(3*n+1)/n * (-2*k/Pgrad)^(1/n)
R^((n+1)/n) = U*(3*n+1)/n * (-2*k/Pgrad)^(1/n)
R = ( U*(3*n+1)/n*(-2*k/Pgrad)^(1/n) )^(n/(n+1))
R = (U*(3*n+1)/n)^(n/(n+1)) * (-2*k/Pgrad)^(1/(n+1))

This is an explicit formula solving for R in terms of your known constants.

When something like this is possible it is usually better to use such a formula rather than methods like bisection to solve it. However, if you are required to do bisection, that is what you will have to do.

[Oliver Jones]

Thanks a lot. Yeah, I had to use the bisection method. Thanks for posting that second comment though, I don't know where I'd gone wrong when rearranging but I was coming out with R=0.6115.... When Matlab was giving me 1.39347.

Once I arranged the formula as you had done, I got the same answer. Used this code by the way if anyone is interested:

function bisection
% Bisection method: Used for solving an equation, and finding an
% approximate solution to an equation.

    clc

    % The number of bisection steps
    n = 32;
    % define the initial interval
    a = input('input a value for lower bound of the interval, a, ');
    b = input('input a value for upper bound of the interval, b, ');
    fprintf('\n initial interval [%g, %g] \n total bisection steps %d\n', a,b,n);

    % Check that a root does exist in the chosen interval
    x_left = a;
    x_right = b;
    f_left = f(x_left);
    f_right = f(x_right);
    if f_left*f_right > 0

    end

    % Bisection: The method
    for i=1:n
        if f_left == 0
            % The exact root of the equation is found at the lower bound
            % of the chosen interval
            fprintf('\n stage %g root %g with zero absolute error \n',i,x_left);
            return;
        end
        if f_right==0
            % The exact root of the equation is found at the upper bound
            % of the chosen interval
            fprintf('\n stage %g root %g with zero absolute error \n',i,x_right);
            return
        end

        % Bisection method: The process
        x_mid = (x_left+x_right)/2.0;
        f_mid = f(x_mid);
        if f_left*f_mid <= 0
            % There is a root in [x_left,x_mid]
            x_right = x_mid;
            f_right = f_mid;
        end
        if f_left*f_mid > 0
            % There is a root in [x_mid,x_right]
            x_left = x_mid;
            f_left = f_mid;
        end

        % Calculate the approximate root for current step
        root = (x_left+x_right)/2.0;
        % Calculate the absolute error for current step
        abs_err=(x_right-x_left)/2.0;
        fprintf('\n stage %g root %g absolute error < %g \n',i,root,abs_err);
    end

    %check satisfaction of equation at end of process
    residual = f(root);
    fprintf('\n final residual = %g \n',residual);
end

    % Subfunction: This defines the equation f(x)= 0
function f_value = f(~)
u = 7;
v = 0;
w = 1;
x = 9;
y = 0;
z = 2;

      k = v + 1;
n = (u./12) + 0.1;
U = 1 + ((w + 1)./((3.*x) + 1));
Pgrad = (-10.*(y + 1))./(z + 1);
f_value = nthroot((U.*((3.*n) + 1))./(n.*(nthroot((Pgrad./(-2.*k)),n))),(1 + (1./n)));
end

[siti bahari]

Hello there

why can't we use the given value of n, U and Pgrad instead of calculating it?

Answer 2

[Walter Roberson]

fun = @(R) -2k * ( ( U * (3n + 1) ) / ( n * R ^ (1 + 1/n) ) ) ^ n - Pgrad

now apply the bisection method to find the R such that fun® = 0

Hint: there is an analytic solution, which is

R = exp((ln(U*(3*n+1)/n)*n-ln(-(1/2)*Pgrad/k))/(n+1))

To as many decimal places as your provide,

R = (20/22776664389) * 7^(6833/16833) * 11^(6833/16833) * 487^(6833/16833) * 4357^(6833/16833) * 2^(15835/16833) * 5^(15835/16833) * 6833^(10000/16833) * 3^(6833/16833) * 239^(6833/16833) * 4649^(6833/16833)