I thought Gaussian quadrature requires that the integral must be transformed to the integral limit of [-1,1]?

Question asked on YouTube: I thought Gaussian quadrature requires that the integral must be transformed to the integral limit of [-1,1]?

The answer is given below.

gaussquadlimits

The document in the above image is given here. This post is brought to you by

Advertisements

An example of Gaussian quadrature rule by using two approaches

Here is an example of using Gaussian quadrature rule through two approaches:

EITHER

by applying it on the original integrand by updating the argument of the integrand

OR

by applying it to the equivalent integrand because of the need to change the limits of integration to: -1 to 1.

http://nm.MathForCollege.com/blog/3pointquadruleexample.pdf

___________________

This post is brought to you by

A MATLAB program to find quadrature points and weights for Gauss-Legendre Quadrature rule

Recently, I got a request how one can find the quadrature and weights of a Gauss-Legendre quadrature rule for large n.  It seems that the internet has these points available free of charge only up to n=12.  Below is the MATLAB program that finds these values for any n.  I tried the program for n=25 and it gave results in a minute or so.  The results output up to 32 significant digits.
_______________________________________________________

% Program to get the quadrature points
% and weight for Gauss-Legendre Quadrature
% Rule
clc
clear all
syms x
% Input n: Quad pt rule
n=14;
% Calculating the Pn(x)
% Legendre Polynomial
% Using recursive relationship
% P(order of polynomial, value of x)
% P(0,x)=1; P(1,x)=0;
% (i+1)*P(i+1,x)=(2*i+1)*x*P(i,x)-i*P(i-1,x)
m=n-1;
P0=1;
P1=x;
for i=1:1:m
    Pn=((2.0*i+1)*x*P1-i*P0)/(i+1.0);
    P0=P1;
    P1=Pn;
end
if n==1
    Pn=P1;
end
Pn=expand(Pn);
quadpts=solve(vpa(Pn,32));
quadpts=sort(quadpts);
% Finding the weights
% Formula for weights is given at
% http://mathworld.wolfram.com/Legendre-GaussQuadrature.html
% Equation (13)
for k=1:1:n
    P0=1;
    P1=x;
    m=n;
    % Calculating P(n+1,x)
    for i=1:1:m
        Pn=((2.0*i+1)*x*P1-i*P0)/(i+1.0);
        P0=P1;
        P1=Pn;
    end
    Pn=P1;
    weights(k)=vpa(2*(1-quadpts(k)^2)/(n+1)^2/ …
                                   subs(Pn,x,quadpts(k))^2,32);
end
    fprintf(‘Quad point rule for n=%g \n’,n)
disp(‘  ‘)
disp(‘Abscissas’)
disp(quadpts)
disp(‘  ‘)
disp(‘Weights’)
disp(weights’)_______________________________________________________ 

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, the textbook on Introduction to Programming Concepts Using MATLAB, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos.  Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.

Using int and solve to find inverse error function in MATLAB

In the previous post, https://autarkaw.wordpress.com/2010/08/24/finding-the-inverse-error-function/, we set up the nonlinear equation to find the inverse of error function.  Using the int and solve MATLAB commands, we write our own program to find the inverse error function.

It is better to download (right click and save target) the program as single quotes in the pasted version do not translate properly when pasted into a mfile editor of MATLAB or you can read the html version for clarity and sample output.

%% FINDING INVERSE ERROR FUNCTION
% In a previous blog at autarkaw.wordpress.com (August 24, 2010),
% we set up a nonlinear equation to find the inverse error function.
% In this blog, we will solve this equation.
% The problem is given at
% http://numericalmethods.eng.usf.edu/blog/inverseerror.pdf
% and we are solving Exercise 1 of the pdf file.

%% TOPIC
% Finding inverse error function

%% SUMMARY

% Language : Matlab 2010a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/inverse_erf_matlab.m;
% Last Revised : August 27, 2010
% Abstract: This program shows you how to find the inverse error function
clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to’)
disp(‘   find the inverse error function’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/inverse_erf_matlab.m’)
disp(‘  ‘)
disp(‘PROBLEM STATEMENT’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/inverseerror.pdf’)
disp(‘        Exercise 1’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   August 27, 2010’)
disp(‘ ‘)

%% INPUTS
% Value of error function
erfx=0.5;

%% DISPLAYING INPUTS

disp(‘INPUTS’)
fprintf(‘ The value of error function= %g’,erfx)
disp(‘  ‘)
disp(‘  ‘)

%% CODE
syms t x
inverse_erf=solve(int(2/sqrt(pi)*exp(-t^2),t,0,x)-erfx);
inverse_erf=double(inverse_erf);
%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘ Value of inverse error function from mfile is= %g’,inverse_erf)
fprintf(‘\n Value of inverse error function using erfinv is= %g’,erfinv(erfx))
disp(‘  ‘)

___________________________________________________

This post is brought to you by

Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu,
the textbook on Numerical Methods with Applications available from the lulu storefront,
the textbook on Introduction to Programming Concepts Using MATLAB, and
the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos

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

Finding the inverse error function

Inverse Error Function

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos and http://www.youtube.com/numericalmethodsguy

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

A video tutorial on Simpson’s 1/3 rule

Simpson’s 1/3 rule is a popular method of conducting numerical integration.  We have recorded a series of short videos on this topic and they are avilable as a playlist at http://www.youtube.com/user/numericalmethodsguy#g/c/2A25C3DC6D8E5616.

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos and http://www.youtube.com/numericalmethodsguy

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

How do I integrate a discrete function in MATLAB?

Many students ask me how do I do this or that in MATLAB.  So I thought why not have a small series of my next few blogs do that.  In this blog, I show you how to integrate a discrete function.

The MATLAB program link is here.

The HTML version of the MATLAB program is here.

_____________________________________________________

%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB.  Most of the questions relate to a mathematical
% procedure.

%% TOPIC
% How do I integrate a discrete function?  Three cases of data are
% discussed.

%% SUMMARY

% Language : MATLAB 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/integrationdiscrete.m;
% Last Revised : April 3, 2009;
% Abstract: This program shows you how to integrate a given discrete function.

clc
clear all

%% INTRODUCTION

disp(‘ABSTRACT’)
disp(‘   This program shows you how to integrate’)
disp(‘   a discrete function’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘   Autar K Kaw of https://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘   http://numericalmethods.eng.usf.edu/blog/integrationdiscrete.m’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘   April 3, 2009’)
disp(‘ ‘)

%% CASE 1

%% INPUTS

% Integrate the discrete function y from x=1 to 6.5
% with y vs x data given as (1,2), (2,7), (4,16), (6.5,18)
% Defining the x-array
x=[1  2  4  6.5];
% Defining the y-array
y=[2  7  16  18];

%% DISPLAYING INPUTS
disp(‘____________________________________’)
disp(‘CASE#1’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION MATCH x(1) AND x(LAST)’)
disp(‘ ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,x(1))
fprintf(‘\n  Upper limit of integration, b= %g’,x(length(x)))
disp(‘ ‘)

%% THE CODE

intvalue=trapz(x,y);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)
disp(‘___________________________________________’)

%% CASE 2

%% INPUTS

% Integrate the discrete function y from x=3 to 6
% with y vs x data given as (1,2), (2,7), (4,16), (6.5,18)
% Defining the x-array
x=[1  2  4  6.5];
% Defining the y-array
y=[2  7  16  18];
% Lower limit of integration, a
a=3;
% Upper limit of integration, b
b=6;
%% DISPLAYING INPUTS

disp(‘CASE#2’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST)’)
disp(‘  ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
% Choose how many divisions you want for splining from a to b
n=1000;
fprintf(‘\n  Number of subdivisions used for splining = %g’,n)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE

xx=a:(b-a)/n:b;
% Using spline to approximate the curve from x(1) to x(last)
yy=spline(x,y,xx);
intvalue=trapz(xx,yy);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)
disp(‘___________________________________________’)
%% CASE 3

%% INPUTS

% Integrate the discrete function y from x=1 to 6.5
% with y vs x data given as (1,2), (4,16), (2,7), (6.5,18)
% The x-data is not in ascending order
% Defining the x-array
x=[1  4   2 6.5];
% Defining the y-array
y=[2  16  7 18];
% Lower limit of integration, a
a=3;
% Upper limit of integration, b
b=6;
%% DISPLAYING INPUTS

disp(‘CASE#3’)
disp(‘LOWER LIMIT AND UPPER LIMITS OF INTEGRATION DO not MATCH x(1) AND x(LAST) ‘)
disp(‘AND X-DATA IS NOT IN ASCENDING OR DESCENDING ORDER’)
disp(‘   ‘)
disp(‘INPUTS’)
disp(‘The x-data is’)
x
disp(‘The y-data is’)
y
fprintf(‘  Lower limit of integration, a= %g’,a)
fprintf(‘\n  Upper limit of integration, b= %g’,b)
% Choose how many divisions you want for splining from a to b
n=1000;
fprintf(‘\n  Number of subdivisions used for splining = %g’,n)
disp(‘  ‘)
disp(‘  ‘)

%% THE CODE
[x,so] = sort(x); % so is the sort order
y = y(so); % y data is now in same order as x data
xx=a:(b-a)/n:b;
% Using spline to approximate the curve from x(1) to x(last)
yy=spline(x,y,xx);
intvalue=trapz(xx,yy);

%% DISPLAYING OUTPUTS

disp(‘OUTPUTS’)
fprintf(‘  Value of integral is = %g’,intvalue)
disp(‘  ‘)

____________________________________________________________

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos and http://www.youtube.com/numericalmethodsguy

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