Numerical integration in Matlab (Simpson's rule)

Question

A cubic polynomial is given by

$$y=\frac{x^3}{6RL}$$

with $R$ and $L$ being constants. Use Matlab and numerical methods to find $x_l$ so that

$$L=\int^{x_l}_0 \sqrt{1+(y')^2} dx$$

when $R=200$ and $L=170.$ Evaluate all the integrals by using the composite Simpson's rule.

Attempt

So if I understand this correctly, we need to use Simpson's rule to evaluate $x_l$ where

$$\int^{x_l}_0 \sqrt{1+(y')^2} dx = \int^{x_l}_0 \sqrt{1+ (\frac{x^2}{2RL})^2} dx= \int^{x_l}_0 \sqrt{1+ (\frac{x^2}{6800})^2} dx = L = 170$$

But I thought the method is used to estimate the integral, not the limits. Here's the code for the composite Simpson's rule:

function I = simpsons(f,a,b,n)
if numel(f)>1 % If the input provided is a vector
    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 % If the input provided is an anonymous function
    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

What do I need to enter as the function? And what do I need to enter for $b$ (the last point of the interval, i.e. $x_l$)?

Any help would be greatly appreciated.

EDIT:

Code for the Bisection method:

a=0; b=170; tol=0.001;
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;
%Bisection method
fa=feval(f,a); fb=feval(f,b);
while abs(b-a)>2*tol
    m=(a+b)/2;
    fm=feval(f,m);
    if fa*fm<=0, b=m;
    else a=m; end
end

Answer 1

For convienience write: $$ I(x)=170-\int_{0}^{x}\sqrt{1+\left(\frac{x^2}{68000}\right)^2}\; dx $$ Put $x_{l}=0$ and $x_{r}=170$. Then $I(x_{l})>0$ and $I(x_{r})<0$, so now you employ the bisection method to find $x_0\in (x_l,x_r)$ such that $I(x_0)=0$ using Simpsons rule to evaluate the integral involved in evaluating $I(x)$ at each step.

Matlab (or rather Gnu-Octave) code for the bisection method:

function [rv]=MyMain(Nknots,tol)
  xx=Bisect(@TT,0,170,Nknots,tol)
end


function [rv]=Simp(ff,a,b,N)
%
%  numerical integration of ff from a to b using
%  Simpson's method with N knots (n should be even)
%
  n=2*ceil(N/2); %working number of knots, forced to even 
  nb2=n/2;
  h=(b-a)/n;
  idxE=2:2:n-1;idxO=1:2:n;

  xE=a+idxE*h;xO=a+idxO*h;
  fE=ff(xE);fO=ff(xO);
  rv=(ff(a)+2*sum(fE)+4*sum(fO)+ff(b))*h/3;
end

function [rv]=FF(x)
  rv=sqrt(1+(x.^2/68000).^2);
end

function [rv]=TT(x,N)
  rv=170-Simp(@FF,0,x,N);
end

function [rv]=Bisect(ff,a,b,Nknots,tol)
%
%  ff is the function that you wand the root of
%  the roots is to be in the interval (a,b) as we
%  require that f(a)*f(b)<0, so a and b bracket a root.
%

  while abs(b-a) > tol
    x = (a+b)/2;
    y = ff(x,Nknots);

    if y == 0 
      rv = x;
      break
    end

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

  rv=(a+b)/2;
end

Of course you could use one of the built in root finders if you are allowed.

Running the computation from the command line:

>> tic,MyMain(100,1e-12),toc
xx =  167.23
Elapsed time is 0.036 seconds.
>>
>> tic,MyMain(10000,1e-12),toc
xx =  167.23
Elapsed time is 0.0530002 seconds.
>>
>> tic,MyMain(10000,1e-8),toc
xx =  167.23
Elapsed time is 0.033006 seconds.
>>
>> tic,MyMain(100,1e-8),toc
xx =  167.23
Elapsed time is 0.0310061 seconds.
>>

Changing from n=100 to n=10000 in Simpson's rule less than doubles the computational time required.

The solution method can be improved somewhat if you send some time analysing the problem. It is quite obvious that the desired root is close to $x=170$, so starting with $0$ as the left hand limit for the bisection can be improved. But really unless you are in a production environment the time spent on further analysis is not really cost effective given how fast this executes as it is.

Answer 2

A couple of notes. Since the integrand is never less than $1$, you know that $x_l<170$. In fact, when $x=170$, the integrand will be about $1.087$, so the eventual answer is bounded below by $\frac{170}{1.087}=156.5$. That gives you a pretty good idea about the starting point.

Bisection is useful when you don't have the derivative available, but here the you have the derivative in the form of the integrand, so it seems to be simpler to implement Newton's method in this case unless bisection is the point of the assignment.

Also you know the integral up to a point near to the actual root at each step after the first, so you need not integrate all the way from zero to the next guess for the root, just from the previous guess. Even if you integrated from zero every time, this program would zip through pretty fast anyway.

EDIT: Newton's method

% simp.m

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

err = 1;
tol = 1.0e-8;
N = 10000;
a = 0;
b = 163;
while abs(err) > tol,
    h = (b-a)/N;
    y = f(a);
    for i = 1:N/2,
        y = y+4*f((2*i-1)*h)+2*f(2*i*h);
    end
    y = y+4*f(b-h)+f(b);
    y = h/3*y-170;
    yp = f(b);
    err = y/yp;
    b = b-err;
end
b