Bisection method and numerical integration

[roam]

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:

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.

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.

 

[Dr Transport]

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...

 

[mfig]

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.

>> 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