Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Bisection method and numerical integration

  1. Apr 14, 2016 #1
    In Matlab I am trying to use the composite Simpson's rule to find ##x_l## so that

    $$170=\int^{x_l}_0 \sqrt{1+(y')^2} dx = \int^{x_l}_0 \sqrt{1+\left( \frac{x^2}{68000} \right)^2} dx $$

    For convenience this can be written as

    $$I(x) = 170 - \int^x_0 \sqrt{1 + (\frac{x^2}{68000})} dx$$

    The limits of integration would be from ##0## to ##x=170##. Now to find the zero of this function I want to employ the bisection method while using Simpsons rule to evaluate the integral involved in evaluating ##I(x)## at each step.

    Here's my code so far:

    Code (Text):

    a=0; b=170;

    for x=[0:]

    f = sqrt(1+((x.^2)./68000).^2);

       %Simpson's rule
       if numel(f)>1
        n=numel(f)-1; h=(b-a)/n;
        I= h/3*(f(1)+2*sum(f(3:2:end-2))+4*sum(f(2:2:end))+f(end));
            else
                h=(b-a)/n; xi=a:h:b;
                I= h/3*(f(xi(1))+2*sum(f(xi(3:2:end-2)))+4*sum(f(xi(2:2:end)))+f(xi(end)));
       end
      f = 170 - I;


    tol=1e-6;
    while abs(b-a) > tol
        x = (a+b)/2;
        y = ff(x);

        if y == 0
            rv = x;
            break
        end

        if ff(a)*y < 0
            b = x;
        else
            a = x;
        end
    end

    rv=(a+b)/2;
    But I get lots of errors and the code does not run. I think it is because I need to get the Simpson's rule to calculate new value for the function at each midpoint of the bisection method. I'm not really sure how to fix this. :confused:

    Alternatively, if bisection does not work with Simpson's method, I would appreciate it if anyone could show me how to exactly incorporate built in root finders into my code.

    Any help would be greatly appreicated.
     
    Last edited: Apr 14, 2016
  2. jcsd
  3. Apr 15, 2016 #2

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    I'd use regula-falsi for the root finding part of this problem. If memory serves me correctly, MATLAB has a simpson's rule built in.....
     
  4. Apr 21, 2016 #3
    The question is confusing because first you say you are trying to find the upper limit of integration ##x_l## (such that the value of the integral is 170), but then later you say that the upper limit of integration is 170.

    If you are trying to find the upper integration limit such that the value of the definite integral is 170, the following will work. quadl is a numerical integrator that comes with MATLAB, but you could also use quadgk or integral to the same effect. Also, fzero is the built-in root finder.

    Code (Text):
    >> f = @(xl) quadl(@(x)sqrt(1+(x.^2/68000).^2),0,xl);
    >> R = fzero(@(x) 170-f(x),4)
    R =
      167.2340
    >> quadl(@(x)sqrt(1+(x.^2/68000).^2),0,R)
    ans =
      170.0000
     
    Last edited: Apr 21, 2016
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Bisection method and numerical integration
  1. Numerical Integration (Replies: 3)

Loading...