Bisection method and numerical integration

  • #1
1,266
11

Main Question or Discussion Point

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

Answers and Replies

  • #2
Dr Transport
Science Advisor
Gold Member
2,329
445
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.....
 
  • #3
275
95
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:
>> 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:

Related Threads on Bisection method and numerical integration

  • Last Post
Replies
3
Views
3K
Replies
1
Views
6K
  • Last Post
Replies
3
Views
26K
Replies
4
Views
61K
Replies
3
Views
9K
Replies
1
Views
2K
Replies
1
Views
2K
  • Last Post
Replies
1
Views
643
  • Last Post
Replies
3
Views
2K
Replies
7
Views
747
Top