Mercurial > forge
changeset 12308:c98bc268ccdc octave-forge
integration: moved to a separate mercurial repo.
| author | carandraug |
|---|---|
| date | Sat, 11 Jan 2014 17:58:39 +0000 |
| parents | 3c7b4702c1f1 |
| children | 3fa7cfe44c4c |
| files | extra/integration/COPYING extra/integration/DESCRIPTION extra/integration/INDEX extra/integration/PKG_ADD extra/integration/PKG_DEL extra/integration/doc/README extra/integration/doc/README.Copying extra/integration/doc/README.gaussq extra/integration/inst/Contents.m extra/integration/inst/count.m extra/integration/inst/cquadnd.m extra/integration/inst/crule.m extra/integration/inst/crule2d.m extra/integration/inst/crule2dgen.m extra/integration/inst/ghrule.m extra/integration/inst/glagrule.m extra/integration/inst/gquad.m extra/integration/inst/gquad2d.m extra/integration/inst/gquad2d6.m extra/integration/inst/gquad2dgen.m extra/integration/inst/gquad6.m extra/integration/inst/gquadnd.m extra/integration/inst/grule.m extra/integration/inst/grule2d.m extra/integration/inst/grule2dgen.m extra/integration/inst/innerfun.m extra/integration/inst/integ1es.m extra/integration/inst/ncrule.m extra/integration/inst/quad2dc.m extra/integration/inst/quad2dcgen.m extra/integration/inst/quad2dg.m extra/integration/inst/quad2dggen.m extra/integration/inst/quadc.m extra/integration/inst/quadg.m extra/integration/inst/quadndg.m extra/integration/inst/test/run.log extra/integration/inst/test/run2dtests.m extra/integration/inst/test/test_ncrule.m extra/integration/inst/test/test_quadg.m extra/integration/inst/test/tests2d.log extra/integration/inst/test/testsnc.log extra/integration/inst/test/testsqg.log extra/integration/inst/testfun/fxpow.m extra/integration/inst/testfun/glimh.m extra/integration/inst/testfun/glimh2.m extra/integration/inst/testfun/gliml.m extra/integration/inst/testfun/gxy.m extra/integration/inst/testfun/gxy1.m extra/integration/inst/testfun/gxy2.m extra/integration/inst/testfun/hx.m extra/integration/inst/testfun/lcrcl.m extra/integration/inst/testfun/lcrcu.m extra/integration/inst/testfun/x25.m extra/integration/inst/testfun/xcubed.m extra/integration/inst/testfun/xsquar.m extra/integration/inst/zero_count.m |
| diffstat | 56 files changed, 0 insertions(+), 2970 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/integration/COPYING Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,337 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 2, June 1991 - - Copyright (C) 1989, 1991 Free Software Foundation, Inc. <http://fsf.org/> - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The licenses for most software are designed to take away your -freedom to share and change it. By contrast, the GNU General Public -License is intended to guarantee your freedom to share and change free -software--to make sure the software is free for all its users. This -General Public License applies to most of the Free Software -Foundation's software and to any other program whose authors commit to -using it. (Some other Free Software Foundation software is covered by -the GNU Library General Public License instead.) You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -this service if you wish), that you receive source code or can get it -if you want it, that you can change the software or use pieces of it -in new free programs; and that you know you can do these things. - - To protect your rights, we need to make restrictions that forbid -anyone to deny you these rights or to ask you to surrender the rights. -These restrictions translate to certain responsibilities for you if you -distribute copies of the software, or if you modify it. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must give the recipients all the rights that -you have. You must make sure that they, too, receive or can get the -source code. And you must show them these terms so they know their -rights. - - We protect your rights with two steps: (1) copyright the software, and -(2) offer you this license which gives you legal permission to copy, -distribute and/or modify the software. - - Also, for each author's protection and ours, we want to make certain -that everyone understands that there is no warranty for this free -software. If the software is modified by someone else and passed on, we -want its recipients to know that what they have is not the original, so -that any problems introduced by others will not reflect on the original -authors' reputations. - - Finally, any free program is threatened constantly by software -patents. We wish to avoid the danger that redistributors of a free -program will individually obtain patent licenses, in effect making the -program proprietary. To prevent this, we have made it clear that any -patent must be licensed for everyone's free use or not licensed at all. - - The precise terms and conditions for copying, distribution and -modification follow. - - GNU GENERAL PUBLIC LICENSE - TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION - - 0. This License applies to any program or other work which contains -a notice placed by the copyright holder saying it may be distributed -under the terms of this General Public License. The "Program", below, -refers to any such program or work, and a "work based on the Program" -means either the Program or any derivative work under copyright law: -that is to say, a work containing the Program or a portion of it, -either verbatim or with modifications and/or translated into another -language. (Hereinafter, translation is included without limitation in -the term "modification".) Each licensee is addressed as "you". - -Activities other than copying, distribution and modification are not -covered by this License; they are outside its scope. The act of -running the Program is not restricted, and the output from the Program -is covered only if its contents constitute a work based on the -Program (independent of having been made by running the Program). -Whether that is true depends on what the Program does. - - 1. You may copy and distribute verbatim copies of the Program's -source code as you receive it, in any medium, provided that you -conspicuously and appropriately publish on each copy an appropriate -copyright notice and disclaimer of warranty; keep intact all the -notices that refer to this License and to the absence of any warranty; -and give any other recipients of the Program a copy of this License -along with the Program. - -You may charge a fee for the physical act of transferring a copy, and -you may at your option offer warranty protection in exchange for a fee. - - 2. You may modify your copy or copies of the Program or any portion -of it, thus forming a work based on the Program, and copy and -distribute such modifications or work under the terms of Section 1 -above, provided that you also meet all of these conditions: - - a) You must cause the modified files to carry prominent notices - stating that you changed the files and the date of any change. - - b) You must cause any work that you distribute or publish, that in - whole or in part contains or is derived from the Program or any - part thereof, to be licensed as a whole at no charge to all third - parties under the terms of this License. - - c) If the modified program normally reads commands interactively - when run, you must cause it, when started running for such - interactive use in the most ordinary way, to print or display an - announcement including an appropriate copyright notice and a - notice that there is no warranty (or else, saying that you provide - a warranty) and that users may redistribute the program under - these conditions, and telling the user how to view a copy of this - License. (Exception: if the Program itself is interactive but - does not normally print such an announcement, your work based on - the Program is not required to print an announcement.) - -These requirements apply to the modified work as a whole. If -identifiable sections of that work are not derived from the Program, -and can be reasonably considered independent and separate works in -themselves, then this License, and its terms, do not apply to those -sections when you distribute them as separate works. But when you -distribute the same sections as part of a whole which is a work based -on the Program, the distribution of the whole must be on the terms of -this License, whose permissions for other licensees extend to the -entire whole, and thus to each and every part regardless of who wrote it. - -Thus, it is not the intent of this section to claim rights or contest -your rights to work written entirely by you; rather, the intent is to -exercise the right to control the distribution of derivative or -collective works based on the Program. - -In addition, mere aggregation of another work not based on the Program -with the Program (or with a work based on the Program) on a volume of -a storage or distribution medium does not bring the other work under -the scope of this License. - - 3. You may copy and distribute the Program (or a work based on it, -under Section 2) in object code or executable form under the terms of -Sections 1 and 2 above provided that you also do one of the following: - - a) Accompany it with the complete corresponding machine-readable - source code, which must be distributed under the terms of Sections - 1 and 2 above on a medium customarily used for software interchange; or, - - b) Accompany it with a written offer, valid for at least three - years, to give any third party, for a charge no more than your - cost of physically performing source distribution, a complete - machine-readable copy of the corresponding source code, to be - distributed under the terms of Sections 1 and 2 above on a medium - customarily used for software interchange; or, - - c) Accompany it with the information you received as to the offer - to distribute corresponding source code. (This alternative is - allowed only for noncommercial distribution and only if you - received the program in object code or executable form with such - an offer, in accord with Subsection b above.) - -The source code for a work means the preferred form of the work for -making modifications to it. For an executable work, complete source -code means all the source code for all modules it contains, plus any -associated interface definition files, plus the scripts used to -control compilation and installation of the executable. However, as a -special exception, the source code distributed need not include -anything that is normally distributed (in either source or binary -form) with the major components (compiler, kernel, and so on) of the -operating system on which the executable runs, unless that component -itself accompanies the executable. - -If distribution of executable or object code is made by offering -access to copy from a designated place, then offering equivalent -access to copy the source code from the same place counts as -distribution of the source code, even though third parties are not -compelled to copy the source along with the object code. - - 4. You may not copy, modify, sublicense, or distribute the Program -except as expressly provided under this License. Any attempt -otherwise to copy, modify, sublicense or distribute the Program is -void, and will automatically terminate your rights under this License. -However, parties who have received copies, or rights, from you under -this License will not have their licenses terminated so long as such -parties remain in full compliance. - - 5. You are not required to accept this License, since you have not -signed it. However, nothing else grants you permission to modify or -distribute the Program or its derivative works. These actions are -prohibited by law if you do not accept this License. Therefore, by -modifying or distributing the Program (or any work based on the -Program), you indicate your acceptance of this License to do so, and -all its terms and conditions for copying, distributing or modifying -the Program or works based on it. - - 6. Each time you redistribute the Program (or any work based on the -Program), the recipient automatically receives a license from the -original licensor to copy, distribute or modify the Program subject to -these terms and conditions. You may not impose any further -restrictions on the recipients' exercise of the rights granted herein. -You are not responsible for enforcing compliance by third parties to -this License. - - 7. If, as a consequence of a court judgment or allegation of patent -infringement or for any other reason (not limited to patent issues), -conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot -distribute so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you -may not distribute the Program at all. For example, if a patent -license would not permit royalty-free redistribution of the Program by -all those who receive copies directly or indirectly through you, then -the only way you could satisfy both it and this License would be to -refrain entirely from distribution of the Program. - -If any portion of this section is held invalid or unenforceable under -any particular circumstance, the balance of the section is intended to -apply and the section as a whole is intended to apply in other -circumstances. - -It is not the purpose of this section to induce you to infringe any -patents or other property right claims or to contest validity of any -such claims; this section has the sole purpose of protecting the -integrity of the free software distribution system, which is -implemented by public license practices. Many people have made -generous contributions to the wide range of software distributed -through that system in reliance on consistent application of that -system; it is up to the author/donor to decide if he or she is willing -to distribute software through any other system and a licensee cannot -impose that choice. - -This section is intended to make thoroughly clear what is believed to -be a consequence of the rest of this License. - - 8. If the distribution and/or use of the Program is restricted in -certain countries either by patents or by copyrighted interfaces, the -original copyright holder who places the Program under this License -may add an explicit geographical distribution limitation excluding -those countries, so that distribution is permitted only in or among -countries not thus excluded. In such case, this License incorporates -the limitation as if written in the body of this License. - - 9. The Free Software Foundation may publish revised and/or new versions -of the General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - -Each version is given a distinguishing version number. If the Program -specifies a version number of this License which applies to it and "any -later version", you have the option of following the terms and conditions -either of that version or of any later version published by the Free -Software Foundation. If the Program does not specify a version number of -this License, you may choose any version ever published by the Free Software -Foundation. - - 10. If you wish to incorporate parts of the Program into other free -programs whose distribution conditions are different, write to the author -to ask for permission. For software which is copyrighted by the Free -Software Foundation, write to the Free Software Foundation; we sometimes -make exceptions for this. Our decision will be guided by the two goals -of preserving the free status of all derivatives of our free software and -of promoting the sharing and reuse of software generally. - - NO WARRANTY - - 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY -FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN -OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES -PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED -OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF -MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS -TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE -PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, -REPAIR OR CORRECTION. - - 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR -REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, -INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING -OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED -TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY -YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER -PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE -POSSIBILITY OF SUCH DAMAGES. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -convey the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - <one line to give the program's name and a brief idea of what it does.> - Copyright (C) <year> <name of author> - - This program is free software; you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation; either version 2 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program; if not, see <http://www.gnu.org/licenses/>. - -Also add information on how to contact you by electronic and paper mail. - -If the program is interactive, make it output a short notice like this -when it starts in an interactive mode: - - Gnomovision version 69, Copyright (C) year name of author - Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. - -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, the commands you use may -be called something other than `show w' and `show c'; they could even be -mouse-clicks or menu items--whatever suits your program. - -You should also get your employer (if you work as a programmer) or your -school, if any, to sign a "copyright disclaimer" for the program, if -necessary. Here is a sample; alter the names: - - Yoyodyne, Inc., hereby disclaims all copyright interest in the program - `Gnomovision' (which makes passes at compilers) written by James Hacker. - - <signature of Ty Coon>, 1 April 1989 - Ty Coon, President of Vice - -This General Public License does not permit incorporating your program into -proprietary programs. If your program is a subroutine library, you may -consider it more useful to permit linking proprietary applications with the -library. If this is what you want to do, use the GNU Library General -Public License instead of this License.
--- a/extra/integration/DESCRIPTION Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,10 +0,0 @@ -Name: Integration -Version: 1.0.7 -Date: 2009-05-03 -Author: Howard Wilson & Bryce Gardner -Maintainer: The Octave Community -Title: Numerical Integration Toolbox -Description: Toolbox for 1-D, 2-D, and n-D Numerical Integration -Depends: octave (>= 2.9.7) -License: GPL version 2 or later -Url: http://octave.sf.net
--- a/extra/integration/INDEX Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -integration >> Integration -Integration - count cquadnd crule crule2d crule2dgen - gquad gquad2d gquad2d6 gquad2dgen gquad6 gquadnd - grule grule2d grule2dgen innerfun ncrule quad2dc - quad2dcgen quad2dg quad2dggen quadc quadg quadndg - zero_count integ1es ghrule glagrule
--- a/extra/integration/PKG_ADD Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,11 +0,0 @@ -# Run this only if the package is installed -if (! exist (fullfile (fileparts (mfilename ("fullpath")), "inst"), "dir")) - dirlist= {"test","testfun"}; - for ii=1:length(dirlist) - addpath (fullfile (fileparts(mfilename("fullpath")),dirlist{ii})) - end -end -# Preserve compatibility with old versions of Octave -if (!exist ("meshgrid","file") && exist ("meshdom","file")) - dispatch("meshgrid","meshdom","any") -end
--- a/extra/integration/PKG_DEL Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -# Run this only if the package is installed -if (! exist (fullfile (fileparts (mfilename ("fullpath")), "inst"), "dir")) - dirlist= {"test","testfun"}; - for ii=1:length(dirlist) - rmpath (fullfile (fileparts(mfilename("fullpath")),dirlist{ii})) - end -end
--- a/extra/integration/doc/README Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,35 +0,0 @@ - Numerical Integration Toolbox - - MATLAB Toolbox for 1-D, 2-D, and n-D Numerical Integration - - The original 1-D routines were obtained from NETLIB and were - written by - Howard Wilson - Department of Engineering Mechanics - University of Alabama - Box 870278 - Tuscaloosa, Alabama 35487-0278 - Phone 205 348-1617 - Email address: HWILSON @ UA1VM.UA.EDU - - The rest of the routines were written by - Bryce Gardner - Ray W. Herrick Laboratories - Purdue University - West Lafayette, IN 47906 - Phone: 317-494-0231 - Fax: 317-494-0787 - Email: gardner@ecn.purdue.edu - - These are the general purpose integration routines: - - quadg.m -- High accuracy replacement for QUAD and QUAD8 (1-D) - quad2dg.m -- 2-D integration over a rectangular region - quad2dggen.m -- 2-D integration over a general region - quadndg.m -- n-D integration over a n-D hyper-rectangular region - - Use the other routines if you want a specific integration quadrature or - specific order of integration quadrature. It is faster if you know that - you need a 10th order integration quadrature to use it directly. If you - use the above routines, the integration will be done with quadratures of - several different orders until the results converge.
--- a/extra/integration/doc/README.Copying Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,59 +0,0 @@ -From help-octave-request@bevo.che.wisc.edu Thu Aug 24 11:45:27 1995 -Received: (from daemon@localhost) by bevo.che.wisc.edu (8.6.12/8.6.12) id LAA14483 for help-octave-outgoing; Thu, 24 Aug 1995 11:45:27 -0500 -Received: from bru.mayo.EDU (bru.mayo.edu [129.176.200.17]) by bevo.che.wisc.edu (8.6.12/8.6.12) with SMTP id LAA14470 for <help-octave@bevo.che.wisc.edu>; Thu, 24 Aug 1995 11:45:26 -0500 -Received: from us0.mayo.EDU by bru.mayo.EDU (4.1/SMI-4.0) - id AA12458; Thu, 24 Aug 95 10:57:03 CDT -Received: from us6.sky2 by us0.mayo.EDU (4.1/SMI-4.1) - id AA08624; Thu, 24 Aug 95 10:57:46 CDT -Date: Thu, 24 Aug 95 10:57:46 CDT -From: vdp@us0.mayo.EDU (Vinayak Dutt) -Received: by us6.sky2 (4.1/SMI-4.1) - id AA18238; Thu, 24 Aug 95 10:58:13 CDT -Subject: Re: Integration toolbox for MATLAB -To: help-octave@bevo.che.wisc.edu -Message-Id: <vdp-9507241558.AA000118199@us6> -Reply-To: Dutt.Vinayak@mayo.EDU -Organization: Ultrasound Research Lab, Mayo Clinic and Foundation -X-Mailer: TkMail-2.0Beta18 -Sender: help-octave-request@bevo.che.wisc.edu - -Hi Octavers: - - I just checked with the author of the original MATLAB toolbox for -gaussian quadrature about using the toolbox with OCTAVE. And its -fine with him if we use/distribute the code for OCTAVE. So the ftp site -maintainers can put the integration toolbox port of mine for -distribution. Here's the message from Bryce Gradner on his permission. - --vinayak -- - ----------- Forwarded message from virgo!autoa.com!bgard@msen.com on Thu, 24 Aug 95 10:15:57 EDT ----------- -Return-Path: <virgo!autoa.com!bgard@msen.com> -Received: from bru.mayo.EDU by us0.mayo.EDU (4.1/SMI-4.1) - id AA08207; Thu, 24 Aug 95 10:37:14 CDT -Received: from heifetz.msen.com by bru.mayo.EDU (4.1/SMI-4.0) - id AA12053; Thu, 24 Aug 95 10:36:28 CDT -Received: from virgo.UUCP by heifetz.msen.com with UUCP - (Smail3.1.28.1 #12) id m0sldyI-0009ZsC; Thu, 24 Aug 95 11:13 EDT -Received: from aquarius.autoa.com by autoa.com (4.1/SMI-4.1) - id AA00684; Thu, 24 Aug 95 10:15:57 EDT -Date: Thu, 24 Aug 95 10:15:57 EDT -From: bgard@autoa.com (Bryce Gardner) -Message-Id: <9508241415.AA00684@autoa.com> -To: Dutt.Vinayak@mayo.EDU -Subject: Re: Integration toolbox for MATLAB - -Hi, - -Sorry it took me so long to get back to you. I am on the octave mailing -list, so I already saw that you were porting the integration routines to -OCTAVE. I am glad to see that you did; it is certainly OK with me to -distribute them with OCTAVE. I'm glad to see that they are being used. - -Bryce - -ps: my email is now bgard@autoa.com - - - -
--- a/extra/integration/doc/README.gaussq Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -% Abstract (March 10,1990) -% -% GAUSSQ - Routines for Composite Gauss Integration -% -% The following routines perform numerical integration by -% composite Gauss integration. A function to give base -% points and weight factors is also included. A document -% comparing results from these formulas and output from the -% function quad provided in MATLAB is provided to illustrate -% the accuracy obtainable using Gauss composite integration. - -% Written by Howard Wilson -% Department of Engineering Mechanics -% University of Alabama -% Box 870278 -% Tuscaloosa, Alabama 35487-0278 -% Phone 205 348-1617 -% Email address: HWILSON @ UA1VM.UA.EDU
--- a/extra/integration/inst/Contents.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,113 +0,0 @@ -% Numerical Integration Toolbox -% -% MATLAB Toolbox for 1-D, 2-D, and n-D Numerical Integration -% -% Edited Version for OCTAVE -% -% The original 1-D routines were obtained from NETLIB and were -% written by -% Howard Wilson -% Department of Engineering Mechanics -% University of Alabama -% Box 870278 -% Tuscaloosa, Alabama 35487-0278 -% Phone 205 348-1617 -% Email address: HWILSON @ UA1VM.UA.EDU -% -% The rest of the routines were written by -% Bryce Gardner -% Ray W. Herrick Laboratories -% Purdue University -% West Lafayette, IN 47906 -% Phone: 317-494-0231 -% Fax: 317-494-0787 -% Email: gardner@ecn.purdue.edu -% -% Easy to use routines: (these routines iteratively integrate with -% higher order quadratures until the integral has -% converged--use these routine unless you want to -% specify the order of integration quadrature that -% is to be used) -% quadg.m -- High accuracy replacement for QUAD and QUAD8 (1-D) -% quad2dg.m -- 2-D integration over a rectangular region -% quad2dggen.m -- 2-D integration over a general region -% quadndg.m -- n-D integration over a n-D hyper-rectangular region -% README.nit -- introductory readme file -% -% The 1-D routines: -% README -- The original readme file by Howard Wilson -% gquad.m -- Integrates a 1-D function with input Gauss -% points and weights (modified by Bryce Gardner to -% handle an optional parameter in the function to be -% integrated and also to calculate the Gauss points -% and weights on the fly) -% gquad6.m -- Integrates a 1-D function with a 6-point quadrature -% grule.m -- Calculates the Gauss points and weights -% run.log -- File with examples -% -% New 1-D routines: -% quadg.m -- High accuracy replacement for QUAD and QUAD8 -% quadc.m -- 1-D Gauss-Chebyshev integration routine -% crule.m -- Calculates the Gauss-Chebyshev points and weights -% ncrule.m -- Calculates the Newton-Coates points and weights -% -% The 2-D routines: -% quad2dg.m -- 2-D integration over a rectangular region -% quad2dc.m -- 2-D integration over a rectangular region with -% a 1/sqrt(1-x.^2)/sqrt(1-y.^2) sinqularity -% gquad2d.m -- Integrates a 2-D function over a square region -% gquad2d6.m -- Integrates a 2-D function over a square region with -% a 6-point quadrature -% quad2dggen.m -- 2-D integration over a general region -% quad2dcgen.m -- 2-D integration over a general region with -% a 1/sqrt(1-x.^2)/sqrt(1-y.^2) sinqularity -% gquad2dgen.m -- Integrates a 2-D function over a variable region -% (That is the limits on the inner integration are -% defined by a function of the variable of integration -% of the outer integral.) -% grule2d.m -- Calculates the Gauss points and weights for gquad2d.m -% grule2dgen.m -- Calculates the Gauss points and weights for -% gquad2dgen.m -% crule2d.m -- Calculates the Gauss-Chebyshev points and weights -% for gquad2d.m -% crule2dgen.m -- Calculates the Gauss-Chebyshev points and weights -% for gquad2dgen.m -% -% The n-D routines: -% quadndg.m -- n-D integration over an n-D hyper-rectangular region -% gquadnd.m -- Integrates a n-D function over -% an n-D hyper-rectangular -% region using a Gauss quadrature -% cquadnd.m -- Integrates a n-D function over -% an n-D hyper-rectangular -% region using a Gauss-Chebyshev quadrature -% innerfun.m -- used internally to gquadnd.m and cquadnd.m -% -% Utility routines: -% count.m -- routine to count the number of function calls -% zero_count.m -- routine to report the number of function calls and -% then to reset the counter -% -% Test scripts: -% run2dtests.m -- 2-D examples and 1-D Gauss-Chebyshev examples -% tests2d.log -- output of run2dtests.m -- Matlab 4.1 on a Sparc 10 -% test_ncrule.m-- m-file to check the Newton-Coates quadrature -% testsnc.log -- output of test_ncrule.m -- Matlab 4.1 on a Sparc 10 -% test_quadg.m -- m-file to check the quadg routine -% testsqg.log -- output of test_quadg.m -- Matlab 4.1 on a Sparc 10 -% -% Test functions: -% xsquar.m -- xsquar(x)=x.^2 -% xcubed.m -- xcubed(x)=x.^3 -% x25.m -- x25(x)=x.^25 -% fxpow.m -- fxpow(x,y)=x.^y -% hx.m -- hx(x)=sum(x.^2) -% gxy.m -- gxy(x,y)=x.^2+y.^2 -% gxy1.m -- gxy1(x,y)=ones(size(x)) -% gxy2.m -- gxy2(x,y)=sqrt(x.^2+y.^2) -% glimh.m -- glimh(y)=3 -% glimh2.m -- glimh(y)=y -% gliml.m -- gliml(y)=0 -% lcrcl.m -- lcrcl(y)=-sqrt(4-y.^2) -% lcrcu.m -- lcrcu(y)=sqrt(4-y.^2) -%
--- a/extra/integration/inst/count.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -global NUM_COUNT - -if ( exist('NUM_COUNT') == 1) - NUM_COUNT=NUM_COUNT+length(x(:)); -else - NUM_COUNT=length(x(:)); -endif
--- a/extra/integration/inst/cquadnd.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,20 +0,0 @@ -function nvol = cquadnd (fun,lowerlim,upperlim,nquad) -%usage: nvol = cquadnd (fun,lowerlim,upperlim,nquad); -% n -- number of dimensions to integrate -% nvol -- value of the n-dimensional integral -% fun -- fun(x) (function to be integrated) in this case treat -% all the different values of x as different variables -% as opposed to different instances of the same variable -% x -- n length vector of coordinates -% lowerlim-- n length vector of lower limits of integration -% upperlim-- n length vector of upper limits of integration -% nquad -- n length vector of number of gauss points -% in each integration - -n=length(lowerlim); -level=n; -x=zeros(n,1); - -nvol = innerfun(fun,lowerlim,upperlim,nquad,n,level,x,'crule'); - -endfunction
--- a/extra/integration/inst/crule.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,32 +0,0 @@ -function [bp,wf]=crule(m) -% -%usage: [bp,wf]=crule(m) -% This function computes Gauss-Chebyshev base points and weight factors -% using the algorithm given by somebody in 'SomeBook', -% page 365, Academic Press, 1975, but modified by a change -% in index variables: j=i+1 and m=n+1. -% The weights are all wf_j=pi/m -% and the base points are bp_j=cos((2j-1)*pi/2/m) -% -% m -- number of Gauss-Chebyshev points (integrates a (2m-1)th order -% polynomial exactly) -% -% The Gauss-Chebyshev Quadrature integrates an integral of the form -% 1 m -% Int ((1/sqrt(1-z^2)) f(z)) dz = pi/m Sum (f(cos((2j-1)*pi/2/m))) -% -1 j=1 -% For compatability with the other Gauss Quadrature routines, I brought -% the weight factor into the summation as -% 1 m -% Int ((1/sqrt(1-z^2)) f(z)) dz = Sum (pi/m * f(cos((2j-1)*pi/2/m))) -% -1 j=1 - -% By Bryce Gardner, Purdue University, Spring 1993. - -j=[1:m]'; - -wf = ones(m,1) * pi / m; - -bp=cos( (2*j-1)*pi / (2*m) ); - -endfunction
--- a/extra/integration/inst/crule2d.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -function [bpx,bpy,wfxy] = crule2d (nquadx,nquady) -% -%usage: [bpx,bpy,wfxy] = crule2d (nquadx,nquady); -% - [bpxv,wfxv]=crule(nquadx); - [bpyv,wfyv]=crule(nquady); - [bpx,bpy]=meshgrid(bpxv,bpyv); - [wfx,wfy]=meshgrid(wfxv,wfyv); -% [bpx,bpy]=meshdom(bpxv,bpyv); -% [wfx,wfy]=meshdom(wfxv,wfyv); - wfxy=wfx.*wfy; -endfunction
--- a/extra/integration/inst/crule2dgen.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function [bpxv,wfxv,bpyv,wfyv] = crule2dgen (nquadx,nquady) -% -%usage: [bpxv,wfxv,bpyv,wfyv] = crule2dgen (nquadx,nquady); -% - [bpxv,wfxv]=crule(nquadx); - [bpyv,wfyv]=crule(nquady); -endfunction
--- a/extra/integration/inst/ghrule.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,67 +0,0 @@ -## Copyright (C) 2003,2011 Eros Sormani, Marco Frontini -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. -## - -## -*- texinfo -*- -## -## @deftypefn {Function File} @ -## {[@var{z}, @var{p}]} = ghrule (@var{n}) -## -## Compute the nodes @var{z} and weights @var{p} -## for the n-point Gauss-Hermite quadrature rule for the -## approximation of the integral of w(x) * f(x) on [-inf,inf] with -## w(x)=exp(-x^2). -## -## Example: -## @example -## [z, p] = ghrule (50); -## abs (sum (p) - quad (@(x) exp (-x.^2), -100, 100, eps)) -## @end example -## -## @seealso{grule,glagrule} -## @end deftypefn - -function [z, p] = ghrule (n) - - if (n <= 1) - z=0; p = sqrt (pi); - return - endif - - jac = zeros (n); - k = [1:n]; - v = sqrt (1/2 * k); - jac += diag (v(1:n-1), 1) + diag (v(1:n-1), -1); - [p, z] = eig (jac); - norm2 =sqrt (diag (p' * p)); % weight normalization - p = (sqrt (pi) * p(1,:)' .^ 2) ./ norm2; % sqrt(pi) = beta0; - z = diag (z); - -endfunction - -%!shared z, p -%! -%!test -%! [z, p] = ghrule (5); -%! err = abs (sum (p) - quad (@(x) exp (-x.^2), -100, 100, eps)); -%! assert (err, 0, sqrt (eps)) -%! -%!test -%! err = abs (dot (z, p)); -%! assert (err, 0, sqrt (eps)) -%! -%!test -%! err = abs (dot (z.^7, p)); -%! assert (err, 0, sqrt (eps)) \ No newline at end of file
--- a/extra/integration/inst/glagrule.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,68 +0,0 @@ -## Copyright (C) 2003,2011 Eros Sormani, Marco Frontini -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; If not, see <http://www.gnu.org/licenses/>. -## - -## -*- texinfo -*- -## -## @deftypefn {Function File} @ -## {[@var{z}, @var{p}]} = glagrule (@var{n}) -## -## Compute the nodes @var{z} and weights @var{p} -## for the n-point Gauss-Laguerre quadrature rule for the -## approximation of the integral of w(x) * f(x) on [0,inf] with -## w(x)=exp(-x). -## -## Example: -## @example -## [z, p] = glagrule (50); -## abs (sum (p) - quad (@(x) exp (-x), 0, 100, eps)) -## @end example -## -## @seealso{grule,ghrule} -## @end deftypefn - -function [z, p] = glagrule (n) - - if (n <= 1) - z=1; - p=1; - return - endif - - jac = zeros (n); - k = [1:n]; - v = sqrt (k.^2); - jac = jac + diag (v(1:n-1), 1) + diag (v(1:n-1), -1) + diag (2*(k-1)+1); - [p,z] = eig (jac); - norm2 = sqrt (diag(p'*p)); % weight normalization - p = (1*p(1,:)'.^2) ./ norm2; % 1 = beta0; - z = diag (z); - -%!shared z, p -%! -%!test -%! [z, p] = glagrule (5); -%! err = abs (sum (p) - quad (@(x) exp (-x), 0, 100, eps)); -%! assert (err, 0, sqrt (eps)) -%! -%!test -%! err = abs (sum (p .* z) - quad (@(x) x .* exp (-x), 0, 100, eps)); -%! assert (err, 0, sqrt (eps)) -%! -%!test -%! err = abs (sum (p .* z.^2) - quad (@(x) x.^2 .* exp (-x), 0, 100, eps)); -%! assert (err, 0, sqrt (eps)) - -
--- a/extra/integration/inst/gquad.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,57 +0,0 @@ -function area = gquad (fun,xlow,xhigh,mparts,bp,wf,y) -% -% area = gquad (fun,xlow,xhigh,mparts,bp,wf) -% or -% area = gquad (fun,xlow,xhigh,mparts,nquad) -% or -% area = gquad (fun,xlow,xhigh,mparts,bp,wf,y) -% This function evaluates the integral of an externally -% defined function fun(x) between limits xlow and xhigh. The -% numerical integration is performed using a composite Gauss -% integration rule. The whole interval is divided into mparts -% subintervals and the integration over each subinterval -% is done with an nquad point Gauss formula which involves base -% points bp and weight factors wf. The normalized interval -% of integration for the bp and wf constants is -1 to +1. The -% algorithm is described by the summation relation -% x=b j=n k=m -% integral( f(x)*dx ) = d1*sum sum( wf(j)*fun(a1+d*k+d1*bp(j)) ) -% x=a j=1 k=1 -% where bp are base points, wf are weight factors -% m = mparts, and n = length(bp) and -% d = (b-a)/m, d1 = d/2, a1 = a-d1 -% The base points and weight factors must first be generated -% by a call to grule of the form [bp,wf] = grule(nquad) -% -% Optional argument, nquad, is used if the Gauss points and weights -% have not been previously calculated. -% -% Optional argument, y, is used if the function, fun is a function -% of x and y. fun(x,y) will be integrated over the range in x for -% the constant, y. - -% by Howard Wilson, U. of Alabama, Spring 1990 -% modified by Bryce Gardner, Purdue U., Spring 1993 to handle -% optional parameter y and also to call with the number of points -% instead of passing the points and weights. - -if ( nargin < 6) - nquad=bp; - [bp,wf]=grule(nquad); -endif -bp=reshape(bp,length(bp),1); wf=reshape(wf,length(wf),1); -d = (xhigh - xlow)/mparts; d2 = d/2; nquad = length(bp); -% x = (d2*bp)*ones([1,mparts]) + (d*ones([nquad,1]))*([1:mparts]); -x1 = (d2.*bp)*ones([1,mparts]) ; -x2 = (d.*ones([nquad,1]))*([1:mparts]); -x = x1 + x2; -x = x + (xlow-d2); -if ( nargin == 7 ) - fv = feval(fun,x,y); -else - fv = feval(fun,x); -endif -wv = wf*ones([1,mparts]); -area=d2.*(sum(wv.*fv)); - -endfunction
--- a/extra/integration/inst/gquad2d.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,63 +0,0 @@ -function vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,b1,b2,w1) -% -%usage: vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy) -% or -% vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,nquadx,nquady) -% This function evaluates the integral of an externally -% defined function fun(x,y) between limits xlow and xhigh -% and ylow and yhigh. The numerical integration is performed -% using a Gauss integration rule. The integration -% is done with an nquadx by nquady Gauss formula which involves base -% point matrices bpx and bpy and weight factor matrix wfxy. The normalized -% interval of integration for the bpx, bpy and wfxy constants is -1 to +1 -% (in x) and -1 to +1 (in y). The algorithm is described by the -% summation relation -% x=b j=nx k=ny -% integral( f(x)*dx ) = J*sum sum( wfxy(j,k)*fun( x(j), y(k) ) ) -% x=a j=1 k=1 -% where wfxy are weight factors, -% nx = nquadx = number of Gauss points in the x-direction, -% ny = nquady = number of Gauss points in the y-direction, -% x = (xhigh-xlow)/2 * bpx + (xhigh+xlow)/2 = mapping function in x, -% y = (yhigh-ylow)/2 * bpy + (yhigh+ylow)/2 = mapping function in y, -% and J = (xhigh-xlow)*(yhigh-ylow)/4 = Jacobian of the mapping. -% The base points and weight factors must first be generated -% by a call to grule of the form [bpx,bpy,wfxy] = grule2d(nquadx,nquady) -% -% The first form of gquad2d is faster when used several times, because -% the points and weights are only calculated once. -% -% The second form of gquad2d is usefull if it is only called once (or a -% few times). - -% by Bryce Gardner, Purdue University, Spring 1993 -% extending Howard Wilson's (U. of Alabama, Spring 1990) -% set of 1-D Gauss quadrature routines. - -if ( nargin == 7 ) - nquadx=b1; - nquady=b2; - [bpx,bpy,wfxy]=grule2d(nquadx,nquady); -elseif ( nargin == 8 ) - bpx = b1; - bpy = b2; - wfxy = w1; -else - disp('Wrong Number of Input Arguments') - return -endif - -%Map to x -qx=(xhigh-xlow)/2; -px=(xhigh+xlow)/2; -x=qx*bpx+px; - -%Map to y -qy=(yhigh-ylow)/2; -py=(yhigh+ylow)/2; -y=qy*bpy+py; - -fv = feval(fun,x,y); -vol = sum(sum(wfxy.*fv))*qx*qy; - -endfunction
--- a/extra/integration/inst/gquad2d6.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,27 +0,0 @@ -function vol = gquad2d6(fun,xlow,xhigh,ylow,yhigh) -% -%usage: vol = gquad6(fun,xlow,xhigh,ylow,yhigh) -% -% ==== Six Point by Six Point Double Integral Gauss Formula ==== -% -% This function determines the volume under an externally -% defined function fun(x,y) between limits xlow and xhigh and -% ylow and yhigh. The numerical integration is performed using -% a gauss integration rule. The integration is done with a -% six point Gauss formula which involves base -% points bpx, bpy and weight factors wfxy. The normalized interval -% of integration for the bp and wf constants is -1 to +1 (in x) and -% -1 to 1 (in y). The algorithm is structured in terms of a -% parameter nquad = 6 which can be changed to accommodate a different -% order formula. - -% by Bryce Gardner, Purdue University, Spring 1993 -% modified from gquad6.m by Howard B. Wilson, U. of Alabama, Spring 1990 - -nquad = 6; -nquadx = nquad; -nquady = nquad; -[bpx,bpy,wfxy] = grule2d(nquadx,nquady); -vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy); - -endfunction
--- a/extra/integration/inst/gquad2dgen.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,64 +0,0 @@ -function vol = gquad2dgen(funxy,limxlow,limxhigh,ylow,yhigh,b1,w1,b2,w2) -% -%usage: vol = gquad2dgen(funxy,limxlow,limxhigh,ylow,yhigh,bpxv,wfxv,bpyv,wfyv); -% or -%usage: vol = gquad2dgen(funxy,limxlow,limxhigh,ylow,yhigh,nquadx,nquady); -% -% This function evaluates a general double integral. The limits of the -% inner integration may be functions of the outer integral's variable of -% integration. Such as -% yhigh ghigh(y) -% Vol = Int Int f(x,y) dx dy -% ylow glow(y) -% where -% funxy = f(x,y) -% limxlow = glow(y) -% limxhigh = ghigh(y) -% and the base points and weighting functions are found from -% [bpxv,wfxv,bpyv,wfyv]=grule2dgen(nquadx,nquady); -% where nquadx and nquady are the number of gauss points in the x- and -% y-directions, respectively. -% -% The first form of gquad2dgen is faster when used several times, because -% the points and weights are only calculated once. -% -% The second form of gquad2dgen is usefull if it is only called once (or a -% few times). - -% by Bryce Gardner, Purdue University, Spring 1993 -% extending Howard Wilson's (U. of Alabama, Spring 1990) -% set of 1-D Gauss quadrature routines to 2-dimensions. - -if ( nargin == 7 ) - nquadx=b1; - nquady=w1; - [bpxv,wfxv,bpyv,wfyv]=grule2dgen(nquadx,nquady); -elseif ( nargin == 9 ) - bpxv = b1; - wfxv = w1; - bpyv = b2; - wfyv = w2; -else - disp('Wrong Number of Input Arguments') - return -endif - -nquady=length(bpyv); -qy=(yhigh-ylow)/2; -py=(yhigh+ylow)/2; - -vol = 0; - -for i=1:nquady - y=qy*bpyv(i)+py; - - xhigh=feval(limxhigh,y); - xlow=feval(limxlow,y); - - vx = gquad(funxy,xlow,xhigh,1,bpxv,wfxv,y); - vol = vol + (wfyv(i) * vx); -endfor - -vol = vol .* qy; - -endfunction
--- a/extra/integration/inst/gquad6.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,43 +0,0 @@ -function area = gquad6(fun,xlow,xhigh,mparts) -% -% -% ==== Six Point Composite Gauss Formula ==== -% ==== With Weight Factors Included ==== -% -% area = gquad6(fun,xlow,xhigh,mparts) -% This function determines the area under an externally -% defined function fun(x) between limits xlow and xhigh. The -% numerical integration is performed using a composite gauss -% integration rule. The whole interval is divided into mparts -% subintervals and the integration over each subinterval -% is done with a six point Gauss formula which involves base -% points bp and weight factors wf. The normalized interval -% of integration for the bp and wf constants is -1 to +1. the -% algorithm is structured in terms of a parameter mquad = 6 which -% can be changed along with bp and wf to accommodate a different -% order formula. The composite algorithm is described by the -% following summation relation -% x=b j=n k=m -% integral( f(x)*dx ) = d1*sum sum( wf(j)*fun(a1+d*k+d1*bp(j)) ) -% x=a j=1 k=1 -% where d = (b-a)/m, d1 = d/2, a1 = a-d1, -% m = mparts, and n = nquad. -% - -% by Howard B. Wilson, U. of Alabama, Spring 1990 - -% The weight factors are -wf = [ 1.71324492379170d-01; 3.60761573048139d-01;... - 4.67913934572691d-01]; wf=[wf;wf([3,2,1])]; - -% The base points are -bp = [-9.32469514203152d-01; -6.61209386466265d-01;... - -2.38619186083197d-01]; bp=[bp;-bp([3,2,1])]; - -d = (xhigh - xlow)/mparts; d2 = d/2; nquad = length(bp); -x = (d2*bp)*ones(1,mparts) + (d*ones(nquad,1))*(1:mparts); -x = x(:) + (xlow-d2); fv=feval(fun,x); wv = wf*ones(1,mparts); - -area=d2*(wv(:)'*fv(:)); - -endfunction
--- a/extra/integration/inst/gquadnd.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,21 +0,0 @@ -function nvol = gquadnd (fun,lowerlim,upperlim,nquad) -% -%usage: nvol = gquadnd (fun,lowerlim,upperlim,nquad); -% n -- number of dimensions to integrate -% nvol -- value of the n-dimensional integral -% fun -- fun(x) (function to be integrated) in this case treat -% all the different values of x as different variables -% as opposed to different instances of the same variable -% x -- n length vector of coordinates -% lowerlim-- n length vector of lower limits of integration -% upperlim-- n length vector of upper limits of integration -% nquad -- n length vector of number of gauss points -% in each integration - -n=length(lowerlim); -level=n; -x=zeros(n,1); - -nvol = innerfun(fun,lowerlim,upperlim,nquad,n,level,x,'grule'); - -endfunction
--- a/extra/integration/inst/grule.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,40 +0,0 @@ -function [bp,wf]=grule(n) -% -% [bp,wf]=grule(n) -% This function computes Gauss base points and weight factors -% using the algorithm given by Davis and Rabinowitz in 'Methods -% of Numerical Integration', page 365, Academic Press, 1975. -% -bp=zeros(n,1); wf=bp; iter=2; m=fix((n+1)/2); e1=n*(n+1); -mm=4*m-1; t=(pi/(4*n+2))*(3:4:mm); nn=(1-(1-1/n)/(8*n*n)); -xo=nn*cos(t); -for j=1:iter - pkm1=1; pk=xo; - for k=2:n - t1=xo.*pk; pkp1=t1-pkm1-(t1-pkm1)/k+t1; - pkm1=pk; pk=pkp1; - endfor - den=1.-xo.*xo; d1=n*(pkm1-xo.*pk); dpn=d1./den; - d2pn=(2.*xo.*dpn-e1.*pk)./den; - d3pn=(4*xo.*d2pn+(2-e1).*dpn)./den; - d4pn=(6*xo.*d3pn+(6-e1).*d2pn)./den; - u=pk./dpn; v=d2pn./dpn; - h=-u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn)))); - p=pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn))); - dp=dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3)); - h=h-p./dp; xo=xo+h; -endfor -bp=-xo-h; -fx=d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(d2pn+(h/4).*(d3pn+(.2*h).*d4pn)))); -wf=2*(1-bp.^2)./(fx.*fx); -if ( (m+m) > n ) - bp(m)=0; -endif -if ( ! ((m+m) == n) ) - m=m-1; -endif -jj=1:m; n1j=(n+1-jj); bp(n1j)=-bp(jj); wf(n1j)=wf(jj); -bp = reshape(bp,length(bp),1); -wf = reshape(wf,length(wf),1); - -endfunction
--- a/extra/integration/inst/grule2d.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -function [bpx,bpy,wfxy] = grule2d (nquadx,nquady) -% -%usage: [bpx,bpy,wfxy] = grule2d (nquadx,nquady); -% - [bpxv,wfxv]=grule(nquadx); - [bpyv,wfyv]=grule(nquady); - [bpx,bpy]=meshgrid(bpxv,bpyv); - [wfx,wfy]=meshgrid(wfxv,wfyv); -% [bpx,bpy]=meshdom(bpxv,bpyv); -% [wfx,wfy]=meshdom(wfxv,wfyv); - wfxy=wfx.*wfy; -endfunction
--- a/extra/integration/inst/grule2dgen.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function [bpxv,wfxv,bpyv,wfyv] = grule2dgen (nquadx,nquady) -% -%usage: [bpxv,wfxv,bpyv,wfyv] = grule2dgen (nquadx,nquady); -% - [bpxv,wfxv]=grule(nquadx); - [bpyv,wfyv]=grule(nquady); -endfunction
--- a/extra/integration/inst/innerfun.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -function int = innerfun(fun,lowerlim,upperlim,nquad,n,level,x,quadrule) -% -%usage: int = innerfun(fun,lowerlim,upperlim,nquad,n,level,x,quadrule); -% - -int = 0; -[bp,wf]=feval(quadrule,nquad(level)); - -xx=x; -qx=(upperlim(level)-lowerlim(level))./2; -px=(upperlim(level)+lowerlim(level))./2; -xlevel=qx.*bp+px; - -nl=nquad(level); -if ( level == 1 ) - for i=1:nl - xx(level)=xlevel(i); - int = int + wf(i) .* feval(fun,xx); - endfor -else - for i=1:nl - xx(level)=xlevel(i); - vint = innerfun(fun,lowerlim,upperlim,nquad,n,level-1,xx,quadrule); - int = int + wf(i) .* vint; - endfor -endif - -int = int .* qx; - -endfunction
--- a/extra/integration/inst/integ1es.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,210 +0,0 @@ -## Copyright (C) Not copyrighted, but placed in the public domain. -## However, it would be appreciated if the statement of authorship -## is preserved. -## -## - -## -*- texinfo -*- -## @deftypefn {Function File} {} integrator(@var{y},@var{x-lower},@var{x-upper}) -## Given a vector of evenly spaced y values on x-lower through -## x-upper, end points included, the numerical integral of y is -## returned. -## -## Unlike "Function" integrators, which need a function routine that -## the integrator can evaluate as it wishes, this is a "Data" -## integrator that enables the integration of the whatever values are -## provided without reference to any other information. Since data -## usually comes in the form of evenly spaced values, this is the -## approach taken. However, this can be used as a "Function" -## integrator as it is in the demo sited below. -## -## This uses a (up through, depending on the number of values -## furnished) ninth order polynomial approximation method and all the -## standard polynomial approximation rules apply hereto. -## -## This is "exact" for 8'th and 9'th orders in x if 9 or more values -## are supplied, 6'th and 7'th for 7 or more, 4'th and 5'th for 5 or -## more, and 2'nd and 3'rd for 3 or more. For higher orders, the -## result will generally become usably accurate with enough data -## points. For illustrations, please see the 'Demo Instructions'. -## -## The function need not be powers of x. Any function that can be -## adequately represented by a ninth order power series will work, -## e.g. exp(-k*x), exp(-k*x^2), sin(k*x), etc. See below. -## -## Demo Instructions: To see how this works with a variety of -## functions, the included demo is run by entering 'demo integ1es', -## without the single quotation marks, at the command line. When the -## source code display, has finished press 'q' to get out of it, and -## then follow the on screen directions. Several test functions are -## available, both with and without experimental error. Note that -## when experimental error, a tolerance of 1%, is included, the error -## in the result, usually in the tenths of a percent, is fairly -## independent of the number of points. -## -## This test serves the purpose, amongst others, of determining how -## many points need be provided to achieve the desired accuracy -## before the data is created. -## -## Beware of trying null value integrals since calculating the -## relative error results in a division by zero. Also, the error for -## e^(ord*x^2) for positive 'ord' is reported as 'NaN' since there is -## no function for comparison. -## @end deftypefn - -## -## Author: Douglas M. Elliott, V1.0, July 1, 2009. -## -## Keywords: numerical integration -## -## - - -function I=integ1es(f,xl,xu); - - - if nargin<3 - disp(" Three inputs are required: the function values, and the lower and upper bounds."); - return - endif - - if xu<xl - disp(" The lower bound must precede the upper."); return - endif - - if xu==xl - disp(" The lower bound can not be the same as the upper."); return - endif - - - C{1}=[1]; - C{2}=[0.5;0.5]; - C{3}=[1;4;1]/3; - C{4}=[3;9;9;3]/8; - C{5}=[14;64;24;64;14]/45; - C{6}=[95;375;250;250;375;95]/288; - C{7}=[41;216;27;272;27;216;41]/140; - C{8}=[751;3577;1323;2989;2989;1323;3577;751]/17280*7; - C{9}=[989;5888;-928;10496;-4540;10496;-928;5888;989]/14175*4; - C{10}=[2857;15741;1080;19344;5778;5778;19344;1080;15741;2857]/89600*9; - C{11}=[16067;106300;-48525;272400;-260550;427368;-260550;272400;-48525; ... - 106300;16067]/299376*5; - C{12}=[2034625;11965622;-1423314;19815805;-4825526;12349588;12349588; ... - -4825526;19815805;-1423314;11965622;2034625]/7257600; - C{13}=[2034625;11965622;-1471442;20354366;-7574721;20812812;-5151324; ... - 20812812;-7574721;20354366;-1471442;11965622;2034625]/7257600; - C{14}=[2034625;11965622;-1471442;20306238;-7036160;18063617;3311900; ... - 3311900;18063617;-7036160;20306238;-1471442;11965622;2034625]/7257600; - C{15}=[2034625;11965622;-1471442;20306238;-7084288;18602178; ... - 562705;11775124;562705;18602178;-7084288;20306238; ... - -1471442;11965622;2034625]/7257600; - C{16}=[2034625;11965622;-1471442;20306238;-7084288;18554050;1101266; ... - 9025929;9025929;1101266;18554050;-7084288;20306238;-1471442; ... - 11965622;2034625]/7257600; - C{17}=[2034625;11965622;-1471442;20306238;-7084288;18554050;1053138; ... - 9564490;6276734;9564490;1053138;18554050;-7084288;20306238; ... - -1471442;11965622;2034625]/7257600; - C{18}=[2034625;11965622;-1471442;20306238;-7084288;18554050;1053138; ... - 9516362;6815295;6815295;9516362;1053138;18554050;-7084288; ... - 20306238;-1471442;11965622;2034625]/7257600; - C{19}=[2034625;11965622;-1471442;20306238;-7084288;18554050; ... - 1053138;9516362;6767167;7353856;6767167;9516362;1053138; ... - 18554050;-7084288;20306238;-1471442;11965622;2034625]/7257600; - C{20}=[2034625;11965622;-1471442;20306238;-7084288;18554050; ... - 1053138;9516362;6767167;7305728;7305728;6767167;9516362; ... - 1053138;18554050;-7084288;20306238;-1471442;11965622;2034625]/7257600; - C10L=[2034625;11965622;-1471442;20306238;-7084288;18554050;1053138; ... - 9516362;6767167;7305728]/7257600; - C10R=[7305728;6767167;9516362;1053138;18554050;-7084288;20306238; ... - -1471442;11965622;2034625]/7257600; - - - g=f; [r,c]=size(g); n=c; if r>c; g=g'; n=r; endif - - if n<21; C10=C{n}; else; C10=[C10L;ones(n-20,1);C10R]; endif - - I=g*C10*(xu-xl)/max(n-1,1); - - -endfunction - - -%!demo -%! -%! clc; input(" Please press 'Enter' to continue. ","s"); -%! % -%! clc; disp(" "); -%! % -%! W$=[" Needs four inputs: the order, the lower and upper"]; -%! W$=[W$,[" bounds, and the integrand function type.\n"]]; -%! X$=[" To add 1% noise, add 0.5 to the integer function number, e.g. 2.5.\n"]; -%! disp(W$); -%! disp("\n The integrand function types are:\n"); -%! disp(" 1. x^ord"); -%! disp(" 2. exp(ord*x)"); -%! disp(" 3. exp(ord*x^2)"); -%! disp(" 4. sin(ord*x)"); -%! disp(" 5. 1/(ord^2 + x^2)\n"); -%! disp(X$); -%! % -%! disp(" "); ord=input(" Enter the order: "); -%! % -%! xl=input(" Enter the lower limit: "); -%! xu=input(" Enter the upper limit: "); -%! if xu<xl -%! disp("\n The lower bound must precede, and can not be greater than the upper.\n"); -%! return -%! endif -%! if xu==xl -%! disp("\n The lower bound can not be the same as the upper.\n"); return -%! endif -%! % -%! fun=input(" Enter the integrand function type: "); fun=abs(round(fun*2))/2; -%! if 5.5<fun -%! disp("\n\n The fourth input, function type, must be 5.5 or less.\n\n"); -%! return -%! endif -%! % -%! % -%! tn=floor(fun); -%! F$={[" x^",num2str(ord)];[" e^(",num2str(ord),"*x)"]; ... -%! [" e^(",num2str(ord),"*x^2)"];[" sin(",num2str(ord),"*x)"]; ... -%! [" 1/(",num2str(ord),"^2+x^2)"];[" x^",num2str(ord),"+1+noise"]}; -%! D$=[" on (",num2str(xl),",",num2str(xu),") "]; -%! if fun~=tn; D$=[", with 1% tolerance,",D$]; endif -%! clc; disp(" "); disp([F${tn},D$]); -%! wn$=" #.##e-014 is just numerical noise.\n"; -%! if fun==tn; disp(wn$); else; disp(" "); endif -%! % -%! % -%! for pts=[2:39,50,100,200,500,1000,10000,100000]; -%! % -%! x=linspace(xl,xu,pts); -%! % -%! if tn==1 -%! f=x.^ord; if fun==1.5; f=f.*(1+(rand(1,pts)-0.5)./50); endif -%! I=integ1es(f,xl,xu); E=I*(ord+1)/(xu^(ord+1)-xl^(ord+1))-1; -%! elseif tn==2 -%! f=exp(ord*x); if fun==2.5; f=f.*(1+(rand(1,pts)-0.5)./50); endif -%! I=integ1es(f,xl,xu); E=I*ord/(exp(xu*ord)-exp(xl*ord))-1; -%! elseif tn==3 -%! f=exp(ord*x.^2); if fun==3.5; f=f.*(1+(rand(1,pts)-0.5)./50); endif -%! I=integ1es(f,xl,xu); -%! if 0<ord; E=NaN; -%! elseif ord==0; E=I/(xu-xl)-1; -%! else; k=sqrt(-ord); E=2*k/pi()^0.5*I/(erf(k*xu)-erf(k*xl))-1; -%! endif -%! elseif tn==4 -%! f=sin(ord*x); if fun==4.5; f=f.*(1+(rand(1,pts)-0.5)./50); endif -%! I=integ1es(f,xl,xu); E=ord*I/(cos(ord*xl)-cos(ord*xu))-1; -%! elseif tn==5 -%! f=1./(ord^2+x.^2); if fun==5.5; f=f.*(1+(rand(1,pts)-0.5)./50); endif -%! I=integ1es(f,xl,xu); E=ord*I/(atan(xu/ord)-atan(xl/ord))-1; -%! endif -%! % -%! printf(" Points = %6.0d I = %12.4e %% RelErr = %9.1e\n",pts,I,E*100); -%! % -%! endfor -%! % -%! % -%! disp(" "); \ No newline at end of file
--- a/extra/integration/inst/ncrule.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,40 +0,0 @@ -function [bp,wf]=ncrule(m) -%usage: [bp,wf]=ncrule(m); -% This function returns the Newton-Coates base points and weight factors -% up to an 8 point Newton-Coates formula. -% -% m -- number of Newton-Coates points (integrates an mth order -% polynomial exactly (or an (m+1)th order for even m)) - -% By Bryce Gardner, Purdue University, Spring 1993. - -if ( m == 1 ) - bp=[-1;1]; - wf=[1;1]; -elseif ( m == 2 ) - bp=[-1;0;1]; - wf=[1;4;1]/3; -elseif ( m == 3 ) - bp=[-1;-1/3;1/3;1]; - wf=[1;3;3;1]/4; -elseif ( m == 4 ) - bp=[-1;-1/2;0;1/2;1]; - wf=[7;32;12;32;7]/45; -elseif ( m == 5 ) - bp=[-1;-3/5;-1/5;1/5;3/5;1]; - wf=[19;75;50;50;75;19]/144; -elseif ( m == 6 ) - bp=[-1;-2/3;-1/3;0;1/3;2/3;1]; - wf=[41;216;27;272;27;216;41]/420; -elseif ( m == 7 ) - bp=[-1;-5/7;-3/7;-1/7;1/7;3/7;5/7;1]; - wf=[751;3577;1323;2989;2989;1323;3577;751]/8640; -else - if ( m != 8 ) - disp('Dont know formula higher than n=8. Returning as if n=8.'); - endif - bp=[-1;-3/4;-1/2;-1/4;0;1/4;1/2;3/4;1]; - wf=[989;5888;-928;10496;-4540;10496;-928;5888;989]/14175; -endif - -endfunction
--- a/extra/integration/inst/quad2dc.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,51 +0,0 @@ -function int = quad2dc(fun,xlow,xhigh,ylow,yhigh,tol) -% -%usage: int = quad2dc('Fun',xlow,xhigh,ylow,yhigh) -%or -% int = quad2dc('Fun',xlow,xhigh,ylow,yhigh,tol) -% -%This function is similar to QUAD or QUAD8 for 2-dimensional integration, -%but it uses a Gaussian-Chebyshev quadrature integration scheme. -% int -- value of the integral -% Fun -- Fun(x,y) (function to be integrated) -% xlow -- lower x limit of integration (should be -xhigh) -% xhigh -- upper x limit of integration -% ylow -- lower y limit of integration (should be -yhigh) -% yhigh -- upper y limit of integration -% tol -- tolerance parameter (optional) -% The Gauss-Chebyshev Quadrature integrates an integral of the form -% yhigh xhigh -% Int ((1/sqrt(1-y^2)) Int ((1/sqrt(1-x^2)) fun(x,y)) dx dy -% -yhigh -xlow - -%This routine could be optimized. - -if ( exist('tol') != 1 ) - tol=1e-3; -elseif ( tol == [] ) - tol=1e-3; -endif - -n=length(xlow); -nquad=2*ones(n,1); -[bpx,bpy,wfxy] = crule2d(2,2); -int_old=gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy); - -converge=0; -for i=1:7 - lim = 2^(i+1); - [bpx,bpy,wfxy] = crule2d(lim,lim); - int=gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy); - - if ( abs(int_old-int) < abs(tol*int) ) - converge=1; - break; - endif - int_old=int; -endfor - -if ( converge == 0 ) - disp('Integral did not converge--singularity likely') -endif - -endfunction
--- a/extra/integration/inst/quad2dcgen.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -function int = quad2dcgen(fun,xlow,xhigh,ylow,yhigh,tol) -% -%usage: int = quad2dcgen('Fun','funxlow','funxhigh',ylow,yhigh) -%or -% int = quad2dcgen('Fun','funxlow','funxhigh',ylow,yhigh,tol) -% -%This function is similar to QUAD or QUAD8 for 2-dimensional integration -%over a general 2-dimensional region, but it uses a Gauss-Chebyshev -%quadrature integration scheme. -%The integral is like: -% yhigh funxhigh(y) -% int = Int (1/sqrt(1-y.^2)) Int (1/sqrt(1-x.^2)) Fun(x,y) dx dy -% ylow funxlow(y) -% -% int -- value of the integral -% Fun -- Fun(x,y) (function to be integrated) -% funxlow -- funxlow(y) -% funxhigh-- funxhigh(y) -% ylow -- lower y limit of integration -% yhigh -- upper y limit of integration -% tol -- tolerance parameter (optional) - -%This routine could be optimized. - -if ( exist('tol') != 1 ) - tol=1e-3; -elseif ( tol==[] ) - tol=1e-3; -endif - -n=length(xlow); -nquad=2*ones(n,1); -[bpxv,wfxv,bpyv,wfyv]=crule2dgen(2,2); -int_old=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,bpxv,wfxv,bpyv,wfyv); - -converge=0; -for i=1:7 - lim = 2^(i+1); - [bpxv,wfxv,bpyv,wfyv]=crule2dgen(lim,lim); - int=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,bpxv,wfxv,bpyv,wfyv); - - if ( abs(int_old-int) < abs(tol*int) ) - converge=1; - break; - endif - int_old=int; -endfor - -if ( converge==0 ) - disp('Integral did not converge--singularity likely') -endif - -endfunction
--- a/extra/integration/inst/quad2dg.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,51 +0,0 @@ -function int = quad2dg(fun,xlow,xhigh,ylow,yhigh,tol) -% -%usage: int = quad2dg('Fun',xlow,xhigh,ylow,yhigh) -%or -% int = quad2dg('Fun',xlow,xhigh,ylow,yhigh,tol) -% -%This function is similar to QUAD or QUAD8 for 2-dimensional integration, -%but it uses a Gaussian quadrature integration scheme. -% int -- value of the integral -% Fun -- Fun(x,y) (function to be integrated) -% xlow -- lower x limit of integration -% xhigh -- upper x limit of integration -% ylow -- lower y limit of integration -% yhigh -- upper y limit of integration -% tol -- tolerance parameter (optional) -%Note that if there are discontinuities the region of integration -%should be broken up into separate pieces. And if there are singularities, -%a more appropriate integration quadrature should be used -%(such as the Gauss-Chebyshev for a specific type of singularity). - -%This routine could be optimized. - -if ( exist('tol') != 1 ) - tol=1e-3; -elseif ( tol==[]) - tol=1e-3; -endif - -n=length(xlow); -nquad=2*ones(n,1); -[bpx,bpy,wfxy] = grule2d(2,2); -int_old=gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy); - -converge=0; -for i=1:7 - lim = 2^(i+1); - [bpx,bpy,wfxy] = grule2d(lim,lim); - int=gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy); - - if ( abs(int_old-int) < abs(tol*int) ) - converge=1; - break; - endif - int_old=int; -endfor - -if ( converge==0 ) - disp('Integral did not converge--singularity likely') -endif - -endfunction
--- a/extra/integration/inst/quad2dggen.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -function int = quad2dggen(fun,xlow,xhigh,ylow,yhigh,tol) -% -%usage: int = quad2dggen('Fun','funxlow','funxhigh',ylow,yhigh) -%or -% int = quad2dggen('Fun','funxlow','funxhigh',ylow,yhigh,tol) -% -%This function is similar to QUAD or QUAD8 for 2-dimensional integration -%over a general 2-dimensional region, but it uses a Gaussian quadrature -%integration scheme. -%The integral is like: -% yhigh funxhigh(y) -% int = Int Int Fun(x,y) dx dy -% ylow funxlow(y) -% -% int -- value of the integral -% Fun -- Fun(x,y) (function to be integrated) -% funxlow -- funxlow(y) -% funxhigh-- funxhigh(y) -% ylow -- lower y limit of integration -% yhigh -- upper y limit of integration -% tol -- tolerance parameter (optional) -%Note that if there are discontinuities the region of integration -%should be broken up into separate pieces. And if there are singularities, -%a more appropriate integration quadrature should be used -%(such as the Gauss-Chebyshev for a specific type of singularity). - -%This routine could be optimized. - -if ( exist('tol') != 1) - tol=1e-3; -elseif ( tol==[] ) - tol=1e-3; -endif - -oldint=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,2,2); - -converge=0; -for i=1:7 - lim = 2^(i+1); - int=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,lim,lim); - - diff = oldint - int; - limit = abs(tol*int); - if ( abs(diff) < limit ) - converge=1; - break; - endif - oldint=int; -endfor - -if ( converge==0 ) - disp('Integral did not converge--singularity likely') -endif - -endfunction
--- a/extra/integration/inst/quadc.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,97 +0,0 @@ -function int = quadc(fun,xlow,xhigh,tol,trace,p1,p2,p3,p4,p5,p6,p7,p8,p9) -% -%usage: int = quadc('Fun',xlow,xhigh) -%or -% int = quadc('Fun',xlow,xhigh,tol) -%or -% int = quadc('Fun',xlow,xhigh,tol,trace,p1,p2,....) -% -%This function works just like QUAD or QUAD8 but uses a Gaussian-Chebyshev -%quadrature integration scheme. -% -% The Gauss-Chebyshev Quadrature integrates an integral of the form -% xhigh -% Int ((1/sqrt(1-x^2)) fun(x)) dx -% -xhigh -%This routine ignores xlow - -global cb2 -global cw2 - -if ( exist('tol') != 1) - tol=1e-3; -elseif ( tol==[] ) - tol=1e-3; -endif -if ( exist('trace') != 1) - trace=0; -elseif ( trace==[] ) - trace=0; -else - trace=1; -endif - -xlow=-xhigh; - -%setup string to call the function -exec_string=['y=',fun,'(x']; -num_parameters=nargin-5; -for i=1:num_parameters - exec_string=[exec_string,',p',int2str(i)]; -endfor -exec_string=[exec_string,');']; - -%setup mapping parameters -jacob=(xhigh-xlow)/2; - -%generate the first two sets of integration points and weights -if ( exist('cb2') != 1 ) - [cb2,cw2]=crule(2); -endif - -x=(cb2+1)*jacob+xlow; -eval(exec_string); -int_old=sum(cw2.*y)*jacob; -if ( trace==1 ) - x_trace=x(:); - y_trace=y(:); -endif - -converge=0; -for i=1:7 - gnum=int2str(2^(i+1)); - vname = ['cb',gnum]; - if ( exist(vname) == 0 ) - estr =['[cb',gnum,',cw',gnum,']=crule(',gnum,');']; - eval(estr); - estr =['global cb' gnum,' cw',gnum,';']; - eval(estr); - endif - estr = ['x=(cb',gnum,'+1)*jacob+xlow;']; - eval(estr); - x=x(:); - eval(exec_string); - estr = ['int=sum(cw',gnum,'.*y)*jacob;']; - eval(estr); - - if ( trace==1 ) - x_trace=[x_trace;x(:)]; - y_trace=[y_trace;y(:)]; - endif - - if ( abs(int_old-int) < abs(tol*int) ) - converge=1; - break; - endif - int_old=int; -endfor - -if ( converge==0 ) - disp('Integral did not converge--singularity likely') -endif - -if ( trace==1 ) - plot(x_trace,y_trace,'+') -endif - -endfunction
--- a/extra/integration/inst/quadg.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,102 +0,0 @@ -function int = quadg(fun,xlow,xhigh,tol,trace,p1,p2,p3,p4,p5,p6,p7,p8,p9) -% -%usage: int = quadg('Fun',xlow,xhigh) -%or -% int = quadg('Fun',xlow,xhigh,tol) -%or -% int = quadg('Fun',xlow,xhigh,tol,trace,p1,p2,....) -% -%This function works just like QUAD or QUAD8 but uses a Gaussian quadrature -%integration scheme. Use this routine instead of QUAD or QUAD8: -% if higher accuracy is desired (this works best if the function, -% 'Fun', can be approximated by a power series) -% or if many similar integrations are going to be done (I think less -% function evaluations will typically be done, but the -% integration points and the weights must be calculated. -% These are saved between integrations so when QUADG -% is called again, the points and weights are all ready -% known.) -% or if the function evaluations are time consuming. -%Note that if there are discontinuities the integral should be broken up into separate -%pieces. And if there are singularities, a more appropriate integration quadrature -%should be used (such as the Gauss-Chebyshev). - -global b2 -global w2 - -if ( exist('tol') != 1 ) - tol=1e-3; -elseif ( tol==[] ) - tol=1e-3; -endif -if ( exist('trace') != 1 ) - trace=0; -elseif ( trace==[] ) - trace=0; -else - trace=1; -endif - -%setup string to call the function -exec_string=['y=',fun,'(x']; -num_parameters=nargin-5; -for i=1:num_parameters - exec_string=[exec_string,',p',int2str(i)]; -endfor -exec_string=[exec_string,');']; - -%setup mapping parameters -jacob=(xhigh-xlow)/2; - -%generate the first two sets of integration points and weights -if ( exist('b2') != 1 ) - [b2,w2]=grule(2); -endif - -x=(b2+1)*jacob+xlow; -eval(exec_string); -int_old=sum(w2.*y)*jacob; -if ( trace==1 ) - x_trace=x(:); - y_trace=y(:); -endif - -converge=0; -for i=1:7 - gnum=int2str(2^(i+1)); - vname = ['b',gnum]; - if ( exist(vname) == 0 ) - estr = ['[b',gnum,',w',gnum,']=grule(',gnum,');']; - eval(estr); - estr = ['global b',gnum,' w',gnum,';']; - eval(estr); - endif - estr = ['x=(b',gnum,'+1)*jacob+xlow;']; - eval(estr); - eval(exec_string); - estr = ['int=sum(w',gnum,'.*y)*jacob;']; - eval(estr); - - if ( trace==1 ) - x_trace=[x_trace;x(:)]; - y_trace=[y_trace;y(:)]; - endif - - if ( abs(int_old-int) < abs(tol*int) ) - converge=1; - break; - endif - int_old=int; -endfor - -if ( converge==0 ) - disp('Integral did not converge--singularity likely') -endif - -if ( trace==1 ) - plot(x_trace,y_trace,'+') -endif - -%gnum,i,length(x_trace) - -endfunction
--- a/extra/integration/inst/quadndg.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,51 +0,0 @@ -function int = quadndg(fun,xlow,xhigh,tol) -% -%usage: int = quadndg('Fun',xlow,xhigh) -%or -% int = quadndg('Fun',xlow,xhigh,tol) -% -%This function is similar to QUAD or QUAD8 for n-dimensional integration, -%but it uses a Gaussian quadrature integration scheme. -% int -- value of the integral -% Fun -- Fun(x) (function to be integrated) in this case treat -% all the different values of x as different variables -% as opposed to different instances of the same variable -% x -- n length vector of coordinates -% xlow -- n length vector of lower limits of integration -% xhigh -- n length vector of upper limits of integration -% tol -- tolerance parameter (optional) -%Note that if there are discontinuities the region of integration -%should be broken up into separate pieces. And if there are singularities, -%a more appropriate integration quadrature should be used -%(such as the Gauss-Chebyshev for a specific type of singularity). - -%This routine could be optimized. - -if ( exist('tol') != 1 ) - tol=1e-3; -elseif ( tol==[] ) - tol=1e-3; -endif - -n=length(xlow); -nquad=2*ones(n,1); -int_old=gquadnd(fun,xlow,xhigh,nquad); - -converge=0; -for i=1:7 - nquad=(2^(i+1))*ones(n,1); - int=gquadnd(fun,xlow,xhigh,nquad); - - if ( abs(int_old-int) < abs(tol*int) ) - converge=1; - break; - endif - int_old=int; -endfor - -if ( converge==0 ) - disp('Integral did not converge--singularity likely') -endif - -endfunction -
--- a/extra/integration/inst/test/run.log Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,90 +0,0 @@ - -% -% === COMPARISON OF RESULTS FROM QUAD, GQUAD, GQUAD6 === - -> format long - -> f=quad('sin',0,pi,1.e-5), % Integrate a simple function - -f = 2.00000006453000 - -> % Use a ten point Gauss formula -> [b10,w10]=grule(10); f=gquad('sin',0,pi,1,b10,w10) - -f = 2.00000000000000 - -> % Try a six point Gauss formula without having to generate -> % base points and weight factors. -> f=gquad6('sin',0,pi,1) - -f = 1.99999999947727, % Accuracy is quite good even with a - % six point formula - -> % Next consider a function having infinite slope at x=0 -> f=3*quad('sqrt',0,1) -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. - -f =1.99998752859620 - -> f=3*gquad('sqrt',0,1,1,b10,w10) - -f = 2.00026812880953, % The exact answer is 2 - -> [b20,w20]=grule(20); f=3*gquad('sqrt',0,1,1,b20,w20) - -f = 2.00003589797091 - -> [b50,w50]=grule(50); f=3*gquad('sqrt',0,1,100,b50,w50) - -f = 2.00000000239878 - -> f=quad('log',0,1) , % Try to integrate a singular function - -Warning: Log of zero - -f = - - -> f= quad('log',eps,1) , % Stop just to the right of x=0 -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. - -f= -1.00421574636997 - -> % How well can a composite Gauss ten point rule do -> % using ten intervals? -> f= gquad('log',0,1,10,b10,w10) - -f= -0.99942637022162, % Exact answer is -1 so results are still - % rather poor. Perhaps a higher order - % will help. - -> [b100,w100]=grule(100); f=gquad('log',0,1,1,b100,w100) - -f= -0.99993748732245, % This result is not much better - -> % Let us try taking a large number of function values -> f=gquad('log',0,1,100,b100,w100) - -f= -0.99999937487322 - -> f=gquad6('log',0,1,1600), % Try the simpler six point formula - -f= -0.99999061950636 - -> % The last several results show that we cannot just 'integrate -> % through' a singularity by using many function values. - -> format short
--- a/extra/integration/inst/test/run2dtests.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,100 +0,0 @@ -format long -% -% x^2+y^2 integrated from -1 to 1 in x and -1 to 1 in y = 8/3 -% -[bx,by,w]=grule2d(2,2); -v=gquad2d('gxy',-1,1,-1,1,bx,by,w) -% or -v=quad2dg('gxy',-1,1,-1,1) -% or -v=gquad2d('gxy',-1,1,-1,1,2,2) -% or -n= gquadnd('hx',[-1;-1],[1;1],[2;2]) -% or -n=quadndg('hx',[-1;-1],[1;1]) -correct_ans=8/3 -% -% x^2+y^2 integrated from 0 to 2 in x and 0 to 2 in y = 32/3 -% -[bx,by,w]=grule2d(2,2); -vol3=gquad2d('gxy',0,2,0,2,bx,by,w) -% or -v=quad2dg('gxy',0,2,0,2) -% or -v=gquad2d('gxy',0,2,0,2,2,2) -% or -n = gquadnd('hx',[0;0],[2;2],[2;2]) -% or -n = quadndg('hx',[0;0],[2;2]) -correct_ans=32/3 -% -% x^2+y^2 integrated from 0 to 3 in x and 0 to 3 in y = 54 -% -[bx,by,w]=grule2d(2,2); -v=gquad2d('gxy',0,3,0,3,bx,by,w) -% or -v=gquad2d('gxy',0,3,0,3,2,2) -% or -n = gquadnd('hx',[0;0],[3;3],[2;2]) -correct_ans=54 -% -% x^2+y^2 integrated from 0 to 3 in x and 0 to 3 in y = 54 -% with the general area of intergration (functional limits in x) -% -[bx2,wx,by2,wy]=grule2dgen(2,2); -v=gquad2dgen('gxy','gliml','glimh',0,3,bx2,wx,by2,wy) -% or -v=gquad2dgen('gxy','gliml','glimh',0,3,2,2) -% or -v=quad2dggen('gxy','gliml','glimh',0,3) -correct_ans=54 -% -% x^2+y^2 integrated from 0 to y in x and 0 to 2 in y = 16/3 -% -v=gquad2dgen('gxy','gliml','glimh2',0,2,bx2,wx,by2,wy) -% or -v=gquad2dgen('gxy','gliml','glimh2',0,2,2,2) -correct_ans=16/3 -% -% 1 integrated from -sqrt(4-y^2) to sqrt(4-y^2) in x -% and -2 to 2 in y = 4*pi -- area of circle with radius 2 or (pi r^2) -% -v=gquad2dgen('gxy1','lcrcl','lcrcu',-2,2,bx2,wx,by2,wy) -% or -v=gquad2dgen('gxy1','lcrcl','lcrcu',-2,2,2,2) -correct_ans=4*pi -% --- same problem better quadratue (more points) -% 1 integrated from -sqrt(4-y^2) to sqrt(4-y^2) in x -% and -2 to 2 in y = 4*pi -- area of circle with radius 2 or (pi r^2) -% -[bx3,wx3,by3,wy3]=grule2dgen(5,5); -v=gquad2dgen('gxy1','lcrcl','lcrcu',-2,2,bx2,wx,by2,wy) -% or -v=gquad2dgen('gxy1','lcrcl','lcrcu',-2,2,5,5) -correct_ans=4*pi -% -% sqrt(x^2+y^2) integrated from -sqrt(4-y^2) to sqrt(4-y^2) in x -% and -2 to 2 in y = 16*pi/3 -% -% Need higher order quadrature -[bx3,wx3,by3,wy3]=grule2dgen(10,10); -v=gquad2dgen('gxy2','lcrcl','lcrcu',-2,2,bx3,wx3,by3,wy3) -% or -v=gquad2dgen('gxy2','lcrcl','lcrcu',-2,2,10,10) -correct_ans=16*pi/3 -% -% 1/sqrt(1-x^2) integrated from -1 to 1 in x = pi -% -% Use Gauss-Chebyshev quadrature -[bpc,wfc]=crule(2); -a=gquad('gxy1',-1,1,1,bpc,wfc) -%or -a=quadc('gxy1',-1,1) -correct_ans=pi -% -% x^2/sqrt(1-x^2) integrated from -1 to 1 in x = pi/2 -% -% Use Gauss-Chebyshev quadrature -a=gquad('xsquar',-1,1,1,bpc,wfc) -a=quadc('xsquar',-1,1) -correct_ans=pi/2
--- a/extra/integration/inst/test/test_ncrule.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -[b1,w1]=ncrule(1); -[b2,w2]=ncrule(2); -[b3,w3]=ncrule(3); -[b4,w4]=ncrule(4); -[b5,w5]=ncrule(5); -[b6,w6]=ncrule(6); -[b7,w7]=ncrule(7); -[b8,w8]=ncrule(8); - -f=gquad('sin',0,pi,1,b1,w1), correct_ans=2 -f=gquad('sin',0,pi,1,b2,w2), correct_ans=2 -f=gquad('sin',0,pi,1,b3,w3), correct_ans=2 -f=gquad('sin',0,pi,1,b4,w4), correct_ans=2 -f=gquad('sin',0,pi,1,b5,w5), correct_ans=2 -f=gquad('sin',0,pi,1,b6,w6), correct_ans=2 -f=gquad('sin',0,pi,1,b7,w7), correct_ans=2 -f=gquad('sin',0,pi,1,b8,w8), correct_ans=2 - -f=gquad('fxpow',0,2,1,b1,w1,1), correct_ans=2 -f=gquad('fxpow',0,2,1,b1,w1,2), correct_ans=8/3 -f=gquad('fxpow',0,2,1,b1,w1,3), correct_ans=4 - -f=gquad('fxpow',0,2,1,b2,w2,2), correct_ans=8/3 -f=gquad('fxpow',0,2,1,b2,w2,3), correct_ans=4 -f=gquad('fxpow',0,2,1,b2,w2,4), correct_ans=32/5 - -f=gquad('fxpow',0,2,1,b3,w3,3), correct_ans=4 -f=gquad('fxpow',0,2,1,b3,w3,4), correct_ans=32/5 -f=gquad('fxpow',0,2,1,b3,w3,5), correct_ans=32/3 - -f=gquad('fxpow',0,2,1,b4,w4,4), correct_ans=32/5 -f=gquad('fxpow',0,2,1,b4,w4,5), correct_ans=32/3 -f=gquad('fxpow',0,2,1,b4,w4,6), correct_ans=128/7 - -f=gquad('fxpow',0,2,1,b5,w5,5),correct_ans=32/3 -f=gquad('fxpow',0,2,1,b5,w5,6),correct_ans=128/7 - -f=gquad('fxpow',0,2,1,b6,w6,6),correct_ans=128/7 -f=gquad('fxpow',0,2,1,b6,w6,7),correct_ans=32 -f=gquad('fxpow',0,2,1,b6,w6,8),correct_ans=512/9 - -f=gquad('fxpow',0,2,1,b7,w7,7),correct_ans=32 -f=gquad('fxpow',0,2,1,b7,w7,8),correct_ans=512/9 - -f=gquad('fxpow',0,2,1,b8,w8,8),correct_ans=512/9 -f=gquad('fxpow',0,2,1,b8,w8,9),correct_ans=512/5 -f=gquad('fxpow',0,2,1,b8,w8,10),correct_ans=2048/11
--- a/extra/integration/inst/test/test_quadg.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,29 +0,0 @@ -f=quad('sin',0,pi), correct_ans=2 -%f=quad8('sin',0,pi), correct_ans=2 -f=quadg('sin',0,pi), correct_ans=2 -f=quadg('sin',0,pi,1e-13), correct_ans=2 % with tighter tolerence - -f=quadg('fxpow',0,2,[],[],1), correct_ans=2 -f=quadg('fxpow',0,2,[],[],2), correct_ans=8/3 - -%zero_count; -%f=quad('fxpow',0,2,[],[],3), correct_ans=4 -%num_functs_in_quad =zero_count -%f=quad8('fxpow',0,2,[],[],3), correct_ans=4 -%num_functs_in_quad8 =zero_count -%f=quadg('fxpow',0,2,[],[],3), correct_ans=4 -%num_functs_in_quadq =zero_count -%disp('quadg has lowest number of function calls') - -%f=quad('fxpow',0,2,[],[],4), correct_ans=32/5 -%f=quad8('fxpow',0,2,[],[],4), correct_ans=32/5 -f=quadg('fxpow',0,2,[],[],4), correct_ans=32/5 - -f=quadg('fxpow',0,2,[],[],5), correct_ans=32/3 -f=quadg('fxpow',0,2,[],[],6), correct_ans=128/7 -f=quadg('fxpow',0,2,[],[],7),correct_ans=32 -f=quadg('fxpow',0,2,[],[],8),correct_ans=512/9 -f=quadg('fxpow',0,2,[],[],9),correct_ans=512/5 - -%f=quad8('fxpow',0,2,[],[],10),correct_ans=2048/11 -f=quadg('fxpow',0,2,[],[],10),correct_ans=2048/11
--- a/extra/integration/inst/test/tests2d.log Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,197 +0,0 @@ ->> run2dtests - -vol1 = - - 2.66666666666667 - - -vol1a = - - 2.66666666666667 - - -vol2 = - - 2.66666666666667 - - -nvol1 = - - 2.66666666666667 - - -nvol1a = - - 2.66666666666667 - - -correct_ans = - - 2.66666666666667 - - -vol3 = - - 10.66666666666667 - - -vol3a = - - 10.66666666666667 - - -vol4 = - - 10.66666666666667 - - -nvol2 = - - 10.66666666666667 - - -nvol2a = - - 10.66666666666667 - - -correct_ans = - - 10.66666666666667 - - -vol5 = - - 54.00000000000002 - - -vol6 = - - 54.00000000000002 - - -correct_ans = - - 54 - - -vol7 = - - 54 - - -vol8 = - - 54 - - -vol8a = - - 54.00000000000001 - - -correct_ans = - - 54 - - -vol9 = - - 5.33333333333333 - - -vol10 = - - 5.33333333333333 - - -correct_ans = - - 5.33333333333333 - - -vol11 = - - 13.06394529484362 - - -vol12 = - - 13.06394529484362 - - -correct_ans = - - 12.56637061435917 - - -vol13 = - - 13.06394529484362 - - -vol14 = - - 12.60725067887481 - - -correct_ans = - - 12.56637061435917 - - -vol15 = - - 16.78129277383199 - - -vol16 = - - 16.78129277383199 - - -correct_ans = - - 16.75516081914556 - - -area1 = - - 3.14159265358979 - - -exec_string = - -y=gxy1(x); - - -area1a = - - 3.14159265358979 - - -correct_ans = - - 3.14159265358979 - - -area2 = - - 1.57079632679490 - - -exec_string = - -y=xsquar(x); - - -area2a = - - 1.57079632679490 - - -correct_ans = - - 1.57079632679490 - ->> diary
--- a/extra/integration/inst/test/testsnc.log Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,302 +0,0 @@ ->> test_ncrule - -f = - - 1.923607162353882e-16 - - -correct_ans = - - 2 - - -f = - - 2.09439510239320 - - -correct_ans = - - 2 - - -f = - - 2.04052428476350 - - -correct_ans = - - 2 - - -f = - - 1.99857073182384 - - -correct_ans = - - 2 - - -f = - - 1.99920309391571 - - -correct_ans = - - 2 - - -f = - - 2.00001781363666 - - -correct_ans = - - 2 - - -f = - - 2.00001086554154 - - -correct_ans = - - 2 - - -f = - - 1.99999983527472 - - -correct_ans = - - 2 - - -f = - - 2 - - -correct_ans = - - 2 - - -f = - - 4 - - -correct_ans = - - 2.66666666666667 - - -f = - - 8 - - -correct_ans = - - 4 - - -f = - - 2.66666666666667 - - -correct_ans = - - 2.66666666666667 - - -f = - - 4 - - -correct_ans = - - 4 - - -f = - - 6.66666666666667 - - -correct_ans = - - 6.40000000000000 - - -f = - - 4 - - -correct_ans = - - 4 - - -f = - - 6.51851851851852 - - -correct_ans = - - 6.40000000000000 - - -f = - - 11.25925925925926 - - -correct_ans = - - 10.66666666666667 - - -f = - - 6.40000000000000 - - -correct_ans = - - 6.40000000000000 - - -f = - - 10.66666666666667 - - -correct_ans = - - 10.66666666666667 - - -f = - - 18.33333333333334 - - -correct_ans = - - 18.28571428571428 - - -f = - - 10.66666666666667 - - -correct_ans = - - 10.66666666666667 - - -f = - - 18.31253333333334 - - -correct_ans = - - 18.28571428571428 - - -f = - - 18.28571428571428 - - -correct_ans = - - 18.28571428571428 - - -f = - - 31.99999999999999 - - -correct_ans = - - 32 - - -f = - - 56.90205761316870 - - -correct_ans = - - 56.88888888888889 - - -f = - - 32.00000000000001 - - -correct_ans = - - 32 - - -f = - - 56.89696413342515 - - -correct_ans = - - 56.88888888888889 - - -f = - - 56.88888888888889 - - -correct_ans = - - 56.88888888888889 - - -f = - - 1.024000000000000e+02 - - -correct_ans = - - 1.024000000000000e+02 - - -f = - - 1.861861979166667e+02 - - -correct_ans = - - 1.861818181818182e+02 - ->> diary
--- a/extra/integration/inst/test/testsqg.log Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,224 +0,0 @@ ->> test_quadg - -f = - - 2.00001659104794 - - -correct_ans = - - 2 - - -f = - - 1.99999999999989 - - -correct_ans = - - 2 - - -f = - - 1.99999999999999 - - -correct_ans = - - 2 - - -f = - - 2.00000000000000 - - -correct_ans = - - 2 - - -f = - - 2 - - -correct_ans = - - 2 - - -f = - - 2.66666666666667 - - -correct_ans = - - 2.66666666666667 - - -f = - - 4 - - -correct_ans = - - 4 - - -num_functs_in_quad = - - 20 - - -f = - - 4 - - -correct_ans = - - 4 - - -num_functs_in_quad8 = - - 33 - - -f = - - 4 - - -correct_ans = - - 4 - - -num_functs_in_quadq = - - 6 - -quadg has lowest number of function calls -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. -Recursion level limit reached in quad. Singularity likely. - -f = - - 6.40003310943535 - - -correct_ans = - - 6.40000000000000 - - -f = - - 6.40000000000000 - - -correct_ans = - - 6.40000000000000 - - -f = - - 6.39999999999999 - - -correct_ans = - - 6.40000000000000 - - -f = - - 10.66666666666666 - - -correct_ans = - - 10.66666666666667 - - -f = - - 18.28571428571427 - - -correct_ans = - - 18.28571428571428 - - -f = - - 31.99999999999996 - - -correct_ans = - - 32 - - -f = - - 56.88888888888881 - - -correct_ans = - - 56.88888888888889 - - -f = - - 1.023999999999999e+02 - - -correct_ans = - - 1.024000000000000e+02 - - -f = - - 1.861818181859950e+02 - - -correct_ans = - - 1.861818181818182e+02 - - -f = - - 1.861818181818181e+02 - - -correct_ans = - - 1.861818181818182e+02 - ->> diary
--- a/extra/integration/inst/testfun/fxpow.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function z=fxpow(x,y) -% -%usage: z=fxpow(x,y) -% - count; - z=x.^y; -endfunction
--- a/extra/integration/inst/testfun/glimh.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function xh=glimh(y) -% -%usage: xh=glimh(y) -% - xh=3; -endfunction
--- a/extra/integration/inst/testfun/glimh2.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function xh=glimh2(y) -% -%usage: xh=glimh2(y) -% - xh=y; -endfunction
--- a/extra/integration/inst/testfun/gliml.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function xl=gliml(y) -% -%usage: xl=gliml(y) -% - xl=0; -endfunction
--- a/extra/integration/inst/testfun/gxy.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,8 +0,0 @@ -function z=gxy(x,y) -% -%usage: z=gxy(x,y); -% -% z=sqrt(x.^2+y.^2); - - z=x.^2+y.^2; -endfunction
--- a/extra/integration/inst/testfun/gxy1.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function z=gxy1(x,y) -% -%usage: z=gxy1(x,y); -% - z=ones(size(x)); - % z=1+x*0+y*0; -endfunction
--- a/extra/integration/inst/testfun/gxy2.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function z=gxy2(x,y) -% -%usage: z=gxy2(x,y); -% - z=sqrt(x.^2+y.^2); - %z=x.^2+y.^2; -endfunction
--- a/extra/integration/inst/testfun/hx.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function z=hx(x) -% -%usage: z=hx(x); -% - z=sum(x.^2); -endfunction
--- a/extra/integration/inst/testfun/lcrcl.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function z=lcrcl(y) -% -%usage: z=lcrcl(y); -% - z=-sqrt(2.^2-y.^2); -endfunction
--- a/extra/integration/inst/testfun/lcrcu.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function z=lcrcu(y) -% -%usage: z=lcrcu(y); -% - z=sqrt(2.^2-y.^2); -endfunction
--- a/extra/integration/inst/testfun/x25.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,7 +0,0 @@ -function y=x25(x) -% -%usage: y=x25(x); -% - y=x.^25; - count -endfunction
--- a/extra/integration/inst/testfun/xcubed.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function y=xcubed(x) -% -%usage: y=xcubed(x) -% - y=x.^3; -endfunction
--- a/extra/integration/inst/testfun/xsquar.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function y=xsquar(x) -% -% x^2 -% - y=x.*x; -endfunction
--- a/extra/integration/inst/zero_count.m Fri Jan 10 14:18:40 2014 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,13 +0,0 @@ -function num_fxpow=zero_count -% -%usage: num_fxpow=zero_count -% -global NUM_COUNT - -if ( exist('NUM_COUNT') == 1 ) - num_fxpow=NUM_COUNT; -endif - -NUM_COUNT=0; - -endfunction