Example to show how numerical ODE solutions can be used to find integrals

In a previous post, I enumerated how we can use numerical ODE techniques like Euler and Runge-Kutta methods to find approximate value of definite integrals. Here is an example. Be sure to do the exercises at the end of the post to appreciate the procedure.

_____________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Advertisements

Can I use numerical solution of ODE techniques to do numerical integration?

Yes.

If you are finding the value of the y=\int_{a}^{b} f(x) dx, then we can solve the integral as an ordinary differential equation as

dy/dx=f(x), y(a)=0

We can now use any of the numerical techniques such as Euler’s methods and Runge-Kutta methods to find the value of y(b) which would be the approximate value of the integral. Use exact techniques of solving linear ODEs with fixed coefficients such as Laplace transforms, and you can have the exact value of the integral.

_____________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Is it just a coincidence – true error in multiple segment Trapezoidal rule gets approximately quartered as the number of segments is doubled?

Look at the table below. This is a table that shows the approximate value of the integral

\int_{8}^{30} 2000 ln \frac{140000}{140000-2100t}-9.8t dt 

as a function of the number of segments used in the Trapezoidal rule and the corresponding true error.

n

Value

Et

1

11868

-807

2

11266

-205

3

11153

-91.4

4

11113

-51.5

5

11094

-33.0

6

11084

-22.9

7

11078

-16.8

8

11074

-12.9

The true error for n=1 is -807 and for n=2 is -205. As you can see the quarter of -807 is approximately -201.75 and close to the true error for n=2. Is this a coincidence?

Look at the true error for n=2 which is -205 and for n=4 is -51.5. As you can see the quarter of -205 is approximately -51.75 and close to the value of the true error for n=4. Is this a coincidence?

No. This is because the true error in a single segment trapezoidal rule is

\frac{(b-a)^3}{12} f^{\prime\prime} (c)

where c is some point not known but in the domain [a,b] of \int_{a}^{b} f(x) dx. It can be then shown (see page 14 of this pdf file for full proof) that for the multiple segment trapezoidal rule, the true error is

\frac{(b-a)^3}{12n^2} f^{\prime\prime}

where the f^{\prime\prime} is an average value of the second derivative of the function f(x) calculated at some point within each of the n segments. Since a and b are constant, and f^{\prime\prime} becomes almost a constant as n increases, the true error is approximately inversely proportional to the square of the number of segments.

Note to the reader: Develop a similar table as given above for an integral of your choice and see it for yourself if the true error gets approximately quartered as the number of segments is doubled.

_____________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Can I use Trapezoidal rule to calculate an improper integral?

For example, \int_{a}^{infinity} e^{-x} dx is an improper integral which can be calculated exactly as e^{-a}.

Can we solve this integral by multiple segment Trapezoidal rule when we already know that the upper limit is infinity?

Yes, we can solve it in spite of the upper limit being infinity. We first need to make a change of variables such as t=\frac{1}{x}, then dx=-\frac{1}{t^2}dt giving \int_{1/a}^{0} e^{-1/t} (-1/t^2) dt

But even for this integral we have issues as the integrand gives division by zero problems at t=0 (although the limit itself is zero at t=0).

The MATLAB program that can be downloaded at http://numericalmethods.eng.usf.edu/blog/trapezoidal_improper.m (better to download it as the single quotes do not translate well with matlab editor) shows that we get accurate results for transformed integral with 256 segments.
% Simulation : Can I use Trapezoidal rule for an improper integral?
% Language : Matlab 2007a
% Authors : Autar Kaw, http://numericalmethods.eng.usf.edu
% mfile available at
% http://numericalmethods.eng.usf.edu/blog/trapezoidal_improper.m
% Last Revised : July 15, 2008
% Abstract: This program shows use of multiple segment Trapezoidal
% rule to integrate exp(-x) from x=a to infinity, where a>0.
clc
clear all
disp(‘This program shows the convergence of getting the value of ‘)
disp(‘an improper integral using multiple segment Trapezoidal rule’)
disp(‘Author: Autar K Kaw. autarkaw.wordpress.com’)

%INPUTS. If you want to experiment, these are the only variables
% you should and can change.
% a = Lower limit of integration
a=1/2;

fprintf(‘\nFinding the integral of exp(-x) from a to infinity where a=%g’,a)

% SIMULATION
syms x
% integrand exp(-x)
f=exp(-x);
valexact=double(int(f,x,a,inf));
fprintf(‘\n\nExact value of integral = %f’,valexact)
disp( ‘ ‘)

% Transformed integrand by change of variables
ft=inline(‘exp(-1/t)*(-1/t^2)’);
% Limits of integration
at=1/a;
b=0;

%finding value of the integral using 16, 32, 64, 128 and 256 segments
for k=4:1:8
n=2^k;
h=(b-at)/n;
sum=0;
for i=1:1:n-1
sum=sum+ft(at+i*h);
end
% See how f(b) is not added as f(b)=infinity. But
% using a different value of the function at one point or
% finite number of points does not change the value of the
% integral. We are substituting f(b)=0
sum=2*sum+ft(at)+0;
sum=(b-at)/(2*n)*sum;
fprintf(‘\nApproximate value of integral =%f with %g segments’,sum,n)
end
disp(‘ ‘)
This program shows the convergence of getting the value of
an improper integral using multiple segment Trapezoidal rule
Author: Autar K Kaw. autarkaw.wordpress.com

Finding the integral of exp(-x) from a to infinity where a=0.5

Exact value of integral = 0.606531

Approximate value of integral =0.605971 with 16 segments
Approximate value of integral =0.606496 with 32 segments
Approximate value of integral =0.606521 with 64 segments
Approximate value of integral =0.606528 with 128 segments
Approximate value of integral =0.606530 with 256 segments

Published with MATLAB® 7.3

Now someone may rightly point out that it works well because the limit of the integrand is zero at the lower limit of integration for the above mentioned integral.

So go ahead, modify the program for doing this improper integral \int_{0}^{4} \frac{1}{\sqrt {x}} dx and see how well it works. In this case the integrand is singular at x=0 and one would have to assume the function to be a number other than infinity at x=0 for applying the multiple segment Trapezoidal rule.

How many segments do you need to get an accurate answer for this integral?

_____________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Experiment for spline interpolation and integration

Background:

The motivation behind the experiment is to understand spline interpolation and numerical integration by finding the volume of water that can be held by a champagne glass.

What does the student do in the lab:

The student chooses one of the odd-shaped champagne glasses (Figure 1). The student measures the outer radius of the champagne glass at different known locations along the height. The student measures the thickness of the glass, so that he/she will be able to find the inner radius of the champagne glass at the locations he/she measured the outer radius. The student pours water to the brim in the champagne glass and checks how much volume the champagne glass holds.

Champagne GlassExercises assigned to the students:
Use MATLAB to solve problems. Use comments, display commands and fprintf statements, sensible variable names and units to explain your work. Staple all the work in the following sequence. Use USCS system of units throughout.

  1. Attach the data sheet on which you collected the data in class.
  2. Find the spline interpolant that curve fits the radius vs height data.
  3. Show the individual points and the spline interpolant of radius vs height on a single plot.
  4. Find how much volume of water the champagne glass would hold.
  5. Compare the above result from problem#4 to the actual volume.
  6. In 100-200 words, type out your conclusions using a word processor. Any formulas should be shown using an equation editor. Any sketches need to be drawn using a drawing software such as Word Drawing. Any plots can be imported from MATLAB.

What materials do you need; where do I buy it; how much do the materials costs?

  1. Champagne Glasses: These glasses, called the Hurricane Plastic Glasses, are available at www.poolsidepineapple.com, part nos. HUR-105, HUR-106, YAR-114. We used glasses made of plastic to avoid breakage. http:/www.poolsidepineapple.com/cart_pages/shopping%20page%20tropical.htm. You can try other places to buy the champagne glasses. About $40 or so for about six pieces including S&H. Better yet, go to a cruise and get souvenir glasses. Whenever you do the experiment, you will remember the good times.
  2. Graduated Cylinder: The graduated cylinder is available at http://scientificsonline.com/, part number 3036286. The cost of the cylinder is $20+S&H.
  3. Vernier Caliper: The caliper is available at http://mcmaster.com part number 20265A49. The cost of the vernier caliper is $60+S&H.
  4. Scale: Need to buy a thin scale for this. Any art-supplies store for $2 or so.

_____________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Length of a curve experiment

In a previous post, I mentioned that I have incorporated experiments in my Numerical Methods course. Here I will discuss the second experiment.

Length of the curve experimentIn this experiment, we find the length of two curves generated from the same points – one curve is a polynomial interpolant and another one is a spline interpolant.

Motivation behind the experiment: In 1901, Runge conducted a numerical experiment to show that higher order interpolation is a bad idea. It was shown that as you use higher order interpolants to approximate f(x)=1/(1+25x2) in [-1,1], the differences between the original function and the interpolants becomes worse. This concept also becomes the basis why we use splines rather than polynomial interpolation to find smooth paths to travel through several discrete points.

What do students do in the lab: A flexible curve (see Figure) of length 12″ made of lead-core construction with graduations in both millimeters and inches is provided. The student needs to draw a curve similar in shape to the Runge’s curve on the provided graphing paper as shown. It just needs to be similar in shape – the student can make the x-domain shorter and the maximum y-value larger or vice-versa. The student just needs to make sure that there is a one-to-one correspondence of values.

Assigned Exercises: Use MATLAB to solve problems (3 thru 6). Use comments, display commands and fprintf statements, sensible variable names and units to explain your work. Staple all the work in the following sequence.

  1. Signed typed affidavit sheet.
  2. Attach the plot you drew in the class. Choose several points (at least nine – do not need to be symmetric) along the curve, including the end points. Write out the co-ordinates on the graphing paper curve as shown in the figure.
  3. Find the polynomial interpolant that curve fits the data. Output the coefficients of the polynomial.
  4. Find the cubic spline interpolant that curve fits the data. Just show the work in the mfile.
  5. Illustrate and show the individual points, polynomial and cubic spline interpolants on a single plot.
  6. Find the length of the two interpolants – the polynomial and the spline interpolant. Calculate the relative difference between the length of each interpolant and the actual length of the flexible curve.
  7. In 100-200 words, type out your conclusions using a word processor. Any formulas should be shown using an equation editor. Any sketches need to be drawn using a drawing software such as Word Drawing. Any plots can be imported from MATLAB.

Where to buy the items for the experiment:

  1. Flexible curves – I bought these via internet at Art City. The brand name is Alvin Tru-Flex Graduated Flexible Curves. Prices range from $5 to $12. Shipping and handling is extra – approximately $6 plus 6% of the price. You may want to buy several 12″ and 16″ flexible curves. I had to send a query to the vendor when I did not receive them within a couple of weeks. Alternatively, call your local Art Store and see if they have them.
  2. Engineering Graph Paper – Staples or Office Depot. Costs about $12 for a pack for 100-sheet pad.
  3. Pencil – Anywhere – My favorite store is the 24-hour Wal-Mart Superstore. $1 for a dozen.
  4. Scale – Anywhere – My favorite store is the 24-hour Wal-Mart Superstore. $1 per unit.

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu

Subscribe to the feed to stay updated and let the information follow you.