[Math] Numerical integration in Matlab (Simpson’s rule)

integrationMATLABnumerical methodssimpsons rule

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

Best Answer

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