From a4bd36217e55c2179f967a40c32eab1d4e520e02 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 27 Jul 2013 17:48:24 +0200 Subject: [PATCH] Integration: simplify namespaces, minor formatting --- src/Examples/Integration.cs | 2 +- src/Numerics/{Integration => }/Integrate.cs | 8 +- .../DoubleExponentialTransformation.cs | 620 ------------------ .../DoubleExponentialTransformation.cs | 620 ++++++++++++++++++ .../NewtonCotesTrapeziumRule.cs | 42 +- .../{Algorithms => }/SimpsonRule.cs | 22 +- src/Numerics/Numerics.csproj | 8 +- .../IntegrationTests/IntegrationTest.cs | 1 - 8 files changed, 661 insertions(+), 662 deletions(-) rename src/Numerics/{Integration => }/Integrate.cs (97%) delete mode 100644 src/Numerics/Integration/Algorithms/DoubleExponentialTransformation.cs create mode 100644 src/Numerics/Integration/DoubleExponentialTransformation.cs rename src/Numerics/Integration/{Algorithms => }/NewtonCotesTrapeziumRule.cs (87%) rename src/Numerics/Integration/{Algorithms => }/SimpsonRule.cs (88%) diff --git a/src/Examples/Integration.cs b/src/Examples/Integration.cs index 99445fa0..43abd834 100644 --- a/src/Examples/Integration.cs +++ b/src/Examples/Integration.cs @@ -25,7 +25,7 @@ // using System; -using MathNet.Numerics.Integration; +using MathNet.Numerics; namespace Examples { diff --git a/src/Numerics/Integration/Integrate.cs b/src/Numerics/Integrate.cs similarity index 97% rename from src/Numerics/Integration/Integrate.cs rename to src/Numerics/Integrate.cs index c3f414b0..5b643adc 100644 --- a/src/Numerics/Integration/Integrate.cs +++ b/src/Numerics/Integrate.cs @@ -28,11 +28,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.Integration -{ - using System; - using Algorithms; +using System; +using MathNet.Numerics.Integration; +namespace MathNet.Numerics +{ /// /// Numeric Integration (Quadrature). /// diff --git a/src/Numerics/Integration/Algorithms/DoubleExponentialTransformation.cs b/src/Numerics/Integration/Algorithms/DoubleExponentialTransformation.cs deleted file mode 100644 index 88aecec0..00000000 --- a/src/Numerics/Integration/Algorithms/DoubleExponentialTransformation.cs +++ /dev/null @@ -1,620 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// http://mathnetnumerics.codeplex.com -// -// Copyright (c) 2009-2010 Math.NET -// -// Permission is hereby granted, free of charge, to any person -// obtaining a copy of this software and associated documentation -// files (the "Software"), to deal in the Software without -// restriction, including without limitation the rights to use, -// copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the -// Software is furnished to do so, subject to the following -// conditions: -// -// The above copyright notice and this permission notice shall be -// included in all copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -// OTHER DEALINGS IN THE SOFTWARE. -// - -namespace MathNet.Numerics.Integration.Algorithms -{ - using System; - using System.Collections.Generic; - - /// - /// Analytic integration algorithm for smooth functions with no discontinuities - /// or derivative discontinuities and no poles inside the interval. - /// - public class DoubleExponentialTransformation - { - #region Precomputed Abcissas and Weights - - /// - /// Precomputed abscissa vector per level. - /// - private static readonly double[][] PrecomputedAbscissas = - new[] - { - new[] - { - 0.00000000000000000000, - 0.95136796407274694574, - 0.99997747719246159286, - 0.99999999999995705839 - }, - new[] - { - 0.67427149224843582608, - 0.99751485645722438683, - 0.99999998887566488198 - }, - new[] - { - 0.37720973816403417379, - 0.85956905868989663517, - 0.98704056050737689169, - 0.99968826402835320905, - 0.99999920473711471266, - 0.99999999995285644818 - }, - new[] - { - 0.19435700332493543161, - 0.53914670538796776905, - 0.78060743898320029925, - 0.91487926326457461091, - 0.97396686819567744856, - 0.99405550663140214329, - 0.99906519645578584642, - 0.99990938469514399984, - 0.99999531604122052843, - 0.99999989278161241838, - 0.99999999914270509218, - 0.99999999999823216531 - }, - new[] - { - 0.097923885287832333262, - 0.28787993274271591456, - 0.46125354393958570440, - 0.61027365750063894488, - 0.73101803479256151149, - 0.82331700550640237006, - 0.88989140278426019808, - 0.93516085752198468323, - 0.96411216422354729193, - 0.98145482667733517003, - 0.99112699244169880223, - 0.99610866543750854254, - 0.99845420876769773751, - 0.99945143443527460584, - 0.99982882207287494166, - 0.99995387100562796075, - 0.99998948201481850361, - 0.99999801714059543208, - 0.99999969889415261122, - 0.99999996423908091534, - 0.99999999678719909830, - 0.99999999978973286224, - 0.99999999999039393352, - 0.99999999999970809734 - }, - new[] - { - 0.049055967305077886315, - 0.14641798429058794053, - 0.24156631953888365838, - 0.33314226457763809244, - 0.41995211127844715849, - 0.50101338937930910152, - 0.57558449063515165995, - 0.64317675898520470128, - 0.70355000514714201566, - 0.75669390863372994941, - 0.80279874134324126576, - 0.84221924635075686382, - 0.87543539763040867837, - 0.90301328151357387064, - 0.92556863406861266645, - 0.94373478605275715685, - 0.95813602271021369012, - 0.96936673289691733517, - 0.97797623518666497298, - 0.98445883116743083087, - 0.98924843109013389601, - 0.99271699719682728538, - 0.99517602615532735426, - 0.99688031812819187372, - 0.99803333631543375402, - 0.99879353429880589929, - 0.99928111192179195541, - 0.99958475035151758732, - 0.99976797159956083506, - 0.99987486504878034648, - 0.99993501992508242369, - 0.99996759306794345976, - 0.99998451990227082442, - 0.99999293787666288565, - 0.99999693244919035751, - 0.99999873547186590954, - 0.99999950700571943689, - 0.99999981889371276701, - 0.99999993755407837378, - 0.99999997987450320175, - 0.99999999396413420165, - 0.99999999832336194826, - 0.99999999957078777261, - 0.99999999989927772326, - 0.99999999997845533741, - 0.99999999999582460688, - 0.99999999999927152627, - 0.99999999999988636130 - }, - new[] - { - 0.024539763574649160379, - 0.073525122985671294475, - 0.12222912220155764235, - 0.17046797238201051811, - 0.21806347346971200463, - 0.26484507658344795046, - 0.31065178055284596083, - 0.35533382516507453330, - 0.39875415046723775644, - 0.44078959903390086627, - 0.48133184611690504422, - 0.52028805069123015958, - 0.55758122826077823080, - 0.59315035359195315880, - 0.62695020805104287950, - 0.65895099174335012438, - 0.68913772506166767176, - 0.71750946748732412721, - 0.74407838354734739913, - 0.76886868676824658459, - 0.79191549237614211447, - 0.81326360850297385168, - 0.83296629391941087564, - 0.85108400798784873261, - 0.86768317577564598669, - 0.88283498824466895513, - 0.89661425428007602579, - 0.90909831816302043511, - 0.92036605303195280235, - 0.93049693799715340631, - 0.93957022393327475539, - 0.94766419061515309734, - 0.95485549580502268541, - 0.96121861515111640753, - 0.96682537031235585284, - 0.97174454156548730892, - 0.97604156025657673933, - 0.97977827580061576265, - 0.98301279148110110558, - 0.98579936302528343597, - 0.98818835380074264243, - 0.99022624046752774694, - 0.99195566300267761562, - 0.99341551316926403900, - 0.99464105571251119672, - 0.99566407681695316965, - 0.99651305464025377317, - 0.99721334704346870224, - 0.99778739195890653083, - 0.99825491617199629344, - 0.99863314864067747762, - 0.99893703483351217373, - 0.99917944893488591716, - 0.99937140114093768690, - 0.99952223765121720422, - 0.99963983134560036519, - 0.99973076151980848263, - 0.99980048143113838630, - 0.99985347277311141171, - 0.99989338654759256426, - 0.99992317012928932869, - 0.99994518061445869309, - 0.99996128480785666613, - 0.99997294642523223656, - 0.99998130127012072679, - 0.99998722128200062811, - 0.99999136844834487344, - 0.99999423962761663478, - 0.99999620334716617675, - 0.99999752962380516793, - 0.99999841381096473542, - 0.99999899541068996962, - 0.99999937270733536947, - 0.99999961398855024275, - 0.99999976602333243312, - 0.99999986037121459941, - 0.99999991800479471056, - 0.99999995264266446185, - 0.99999997311323594362, - 0.99999998500307631173, - 0.99999999178645609907, - 0.99999999558563361584, - 0.99999999767323673790, - 0.99999999879798350040, - 0.99999999939177687583, - 0.99999999969875436925, - 0.99999999985405611550, - 0.99999999993088839501, - 0.99999999996803321674, - 0.99999999998556879008, - 0.99999999999364632387, - 0.99999999999727404948, - 0.99999999999886126543, - 0.99999999999953722654, - 0.99999999999981720098, - 0.99999999999992987953 - } - }; - - /// - /// Precomputed weight vector per level. - /// - private static readonly double[][] PrecomputedWeights = - new[] - { - new[] - { - 1.5707963267948966192, - 0.230022394514788685, - 0.00026620051375271690866, - 1.3581784274539090834e-12 - }, - new[] - { - 0.96597657941230114801, - 0.018343166989927842087, - 2.1431204556943039358e-7 - }, - new[] - { - 1.3896147592472563229, - 0.53107827542805397476, - 0.076385743570832304188, - 0.0029025177479013135936, - 0.000011983701363170720047, - 1.1631165814255782766e-9 - }, - new[] - { - 1.5232837186347052132, - 1.1934630258491569639, - 0.73743784836154784136, - 0.36046141846934367417, - 0.13742210773316772341, - 0.039175005493600779072, - 0.0077426010260642407123, - 0.00094994680428346871691, - 0.000062482559240744082891, - 1.8263320593710659699e-6, - 1.8687282268736410132e-8, - 4.9378538776631926964e-11 - }, - new[] - { - 1.5587733555333301451, - 1.466014426716965781, - 1.297475750424977998, - 1.0816349854900704074, - 0.85017285645662006895, - 0.63040513516474369106, - 0.44083323627385823707, - 0.290240679312454185, - 0.17932441211072829296, - 0.10343215422333290062, - 0.055289683742240583845, - 0.027133510013712003219, - 0.012083543599157953493, - 0.0048162981439284630173, - 0.0016908739981426396472, - 0.00051339382406790336017, - 0.00013205234125609974879, - 0.000028110164327940134749, - 4.8237182032615502124e-6, - 6.4777566035929719908e-7, - 6.5835185127183396672e-8, - 4.8760060974240625869e-9, - 2.5216347918530148572e-10, - 8.6759314149796046502e-12 - }, - new[] - { - 1.5677814313072218572, - 1.5438811161769592204, - 1.4972262225410362896, - 1.4300083548722996676, - 1.3452788847662516615, - 1.2467012074518577048, - 1.1382722433763053734, - 1.0240449331118114483, - 0.90787937915489531693, - 0.79324270082051671787, - 0.68306851634426375464, - 0.57967810308778764708, - 0.48475809121475539287, - 0.39938474152571713515, - 0.32408253961152890402, - 0.258904639514053516, - 0.20352399885860174519, - 0.15732620348436615027, - 0.11949741128869592428, - 0.089103139240941462841, - 0.065155533432536205042, - 0.046668208054846613644, - 0.032698732726609031113, - 0.022379471063648476483, - 0.014937835096050129696, - 0.0097072237393916892692, - 0.0061300376320830301252, - 0.0037542509774318343023, - 0.0022250827064786427022, - 0.0012733279447082382027, - 0.0007018595156842422708, - 0.00037166693621677760301, - 0.00018856442976700318572, - 0.000091390817490710122732, - 0.000042183183841757600604, - 0.000018481813599879217116, - 7.6595758525203162562e-6, - 2.9916615878138787094e-6, - 1.0968835125901264732e-6, - 3.7595411862360630091e-7, - 1.1992442782902770219e-7, - 3.5434777171421953043e-8, - 9.6498888961089633609e-9, - 2.4091773256475940779e-9, - 5.482835779709497755e-10, - 1.1306055347494680536e-10, - 2.0989335404511469109e-11, - 3.4841937670261059685e-12 - }, - new[] - { - 1.5700420292795931467, - 1.5640214037732320999, - 1.5520531698454121192, - 1.5342817381543034316, - 1.5109197230741697127, - 1.48224329788553807, - 1.4485862549613225916, - 1.4103329714462590129, - 1.3679105116808964881, - 1.3217801174437728579, - 1.2724283455378627082, - 1.2203581095793582207, - 1.1660798699324345766, - 1.1101031939653403796, - 1.0529288799552666556, - 0.99504180404613271514, - 0.93690461274566793366, - 0.87895234555278212039, - 0.82158803526696470334, - 0.7651792989089561367, - 0.71005590120546898385, - 0.65650824613162753076, - 0.60478673057840362158, - 0.55510187800363350959, - 0.5076251588319080997, - 0.4624903980553677613, - 0.41979566844501548066, - 0.37960556938665160999, - 0.3419537959230168323, - 0.30684590941791694932, - 0.27426222968906810637, - 0.24416077786983990868, - 0.21648020911729617038, - 0.19114268413342749532, - 0.16805663794826916233, - 0.14711941325785693248, - 0.12821973363120098675, - 0.11123999898874453035, - 0.096058391865189467849, - 0.082550788110701737654, - 0.070592469906866999352, - 0.060059642358636300319, - 0.05083075757257047107, - 0.042787652157725676034, - 0.035816505604196436523, - 0.029808628117310126969, - 0.024661087314753282511, - 0.020277183817500123926, - 0.016566786254247575375, - 0.013446536605285730674, - 0.010839937168255907211, - 0.0086773307495391815854, - 0.0068957859690660035329, - 0.0054388997976239984331, - 0.0042565295990178580165, - 0.0033044669940348302363, - 0.0025440657675291729678, - 0.0019418357759843675814, - 0.0014690143599429791058, - 0.0011011261134519383862, - 0.00081754101332469493115, - 0.00060103987991147422573, - 0.00043739495615911687786, - 0.00031497209186021200274, - 0.00022435965205008549104, - 0.00015802788400701191949, - 0.00011002112846666697224, - 0.000075683996586201477788, - 0.000051421497447658802092, - 0.0000344921247593431977, - 0.000022832118109036146591, - 0.000014908514031870608449, - 9.5981941283784710776e-6, - 6.0899100320949039256e-6, - 3.8061983264644899045e-6, - 2.3421667208528096843e-6, - 1.4183067155493917523e-6, - 8.4473756384859863469e-7, - 4.9458288702754198508e-7, - 2.8449923659159806339e-7, - 1.6069394579076224911e-7, - 8.9071395140242387124e-8, - 4.8420950198072369669e-8, - 2.579956822953589238e-8, - 1.3464645522302038796e-8, - 6.8784610955899001111e-9, - 3.4371856744650090511e-9, - 1.6788897682161906807e-9, - 8.0099784479729665356e-10, - 3.7299501843052790038e-10, - 1.6939457789411646876e-10, - 7.4967397573818224522e-11, - 3.230446433325236576e-11, - 1.3542512912336274432e-11, - 5.5182369468174885821e-12, - 2.1835922099233609052e-12 - } - }; - - #endregion - - /// - /// Maximum number of iterations, until the asked - /// maximum error is (likely to be) satisfied. - /// - private const int NumberOfMaximumLevels = 10; - - /// - /// Abscissa vector per level provider. - /// - private IEnumerable _levelAbcissas; - - /// - /// Weight vector per level provider. - /// - private IEnumerable _levelWeights; - - /// - /// Approximate the integral by the double exponential transformation - /// - /// The analytic smooth function to integrate. - /// Where the interval starts, inclusive and finite. - /// Where the interval stops, inclusive and finite. - /// The expected relative accuracy of the approximation. - /// Approximation of the finite integral in the given interval. - public double Integrate( - Func f, - double intervalBegin, - double intervalEnd, - double targetRelativeError) - { - if (_levelAbcissas == null) - { - _levelAbcissas = ProvideLevelAbcissas(); - _levelWeights = ProvideLevelWeights(); - } - - return NewtonCotesTrapeziumRule.IntegrateAdaptiveTransformedOdd( - f, - intervalBegin, - intervalEnd, - _levelAbcissas, - _levelWeights, - 1, - targetRelativeError); - } - - /// - /// Abscissa vector per level provider. - /// - /// Level Enumerator. - internal static IEnumerable ProvideLevelAbcissas() - { - for (int i = 0; i < NumberOfMaximumLevels; i++) - { - yield return EvaluateAbcissas(i); - } - } - - /// - /// Weight vector per level provider. - /// - /// Level Enumerator. - internal static IEnumerable ProvideLevelWeights() - { - for (int i = 0; i < NumberOfMaximumLevels; i++) - { - yield return EvaluateWeights(i); - } - } - - /// - /// Compute the abscissa vector for a single level. - /// - /// The level to evaluate the abscissa vector for. - /// Abscissa Vector. - private static double[] EvaluateAbcissas(int level) - { - if (level < PrecomputedAbscissas.Length) - { - return PrecomputedAbscissas[level]; - } - - double step = level <= 1 ? 1.0 : (1.0 / (2 << (level - 2))); - double offset = level == 0 ? 0.0 : (1.0 / (2 << (level - 1))); - int length = level == 0 ? 4 : (3 << (level - 1)); - - double t = 0; - double[] abcissas = new double[length]; - for (int i = 0; i < abcissas.Length; i++) - { - double arg = offset + t; - t += step; - - abcissas[i] = Math.Tanh(Constants.PiOver2 * Math.Sinh(arg)); - } - - return abcissas; - } - - /// - /// Compute the weight vector for a single level. - /// - /// The level to evaluate the weight vector for. - /// Weight Vector. - private static double[] EvaluateWeights(int level) - { - if (level < PrecomputedWeights.Length) - { - return PrecomputedWeights[level]; - } - - double step = level <= 1 ? 1.0 : (1.0 / (2 << (level - 2))); - double offset = level == 0 ? 0.0 : (1.0 / (2 << (level - 1))); - int length = level == 0 ? 4 : (3 << (level - 1)); - - double t = 0; - double[] weights = new double[length]; - for (int i = 0; i < weights.Length; i++) - { - double arg = offset + t; - t += step; - - // TODO: reuse abcissas as computed in EvaluateAbcissas - double abcissa = Math.Tanh(Constants.PiOver2 * Math.Sinh(arg)); - weights[i] = Constants.PiOver2 * (1 - (abcissa * abcissa)) * Math.Cosh(arg); - } - - return weights; - } - } -} \ No newline at end of file diff --git a/src/Numerics/Integration/DoubleExponentialTransformation.cs b/src/Numerics/Integration/DoubleExponentialTransformation.cs new file mode 100644 index 00000000..da46e769 --- /dev/null +++ b/src/Numerics/Integration/DoubleExponentialTransformation.cs @@ -0,0 +1,620 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2010 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; + +namespace MathNet.Numerics.Integration +{ + /// + /// Analytic integration algorithm for smooth functions with no discontinuities + /// or derivative discontinuities and no poles inside the interval. + /// + public class DoubleExponentialTransformation + { + #region Precomputed Abcissas and Weights + + /// + /// Precomputed abscissa vector per level. + /// + static readonly double[][] PrecomputedAbscissas = + new[] + { + new[] + { + 0.00000000000000000000, + 0.95136796407274694574, + 0.99997747719246159286, + 0.99999999999995705839 + }, + new[] + { + 0.67427149224843582608, + 0.99751485645722438683, + 0.99999998887566488198 + }, + new[] + { + 0.37720973816403417379, + 0.85956905868989663517, + 0.98704056050737689169, + 0.99968826402835320905, + 0.99999920473711471266, + 0.99999999995285644818 + }, + new[] + { + 0.19435700332493543161, + 0.53914670538796776905, + 0.78060743898320029925, + 0.91487926326457461091, + 0.97396686819567744856, + 0.99405550663140214329, + 0.99906519645578584642, + 0.99990938469514399984, + 0.99999531604122052843, + 0.99999989278161241838, + 0.99999999914270509218, + 0.99999999999823216531 + }, + new[] + { + 0.097923885287832333262, + 0.28787993274271591456, + 0.46125354393958570440, + 0.61027365750063894488, + 0.73101803479256151149, + 0.82331700550640237006, + 0.88989140278426019808, + 0.93516085752198468323, + 0.96411216422354729193, + 0.98145482667733517003, + 0.99112699244169880223, + 0.99610866543750854254, + 0.99845420876769773751, + 0.99945143443527460584, + 0.99982882207287494166, + 0.99995387100562796075, + 0.99998948201481850361, + 0.99999801714059543208, + 0.99999969889415261122, + 0.99999996423908091534, + 0.99999999678719909830, + 0.99999999978973286224, + 0.99999999999039393352, + 0.99999999999970809734 + }, + new[] + { + 0.049055967305077886315, + 0.14641798429058794053, + 0.24156631953888365838, + 0.33314226457763809244, + 0.41995211127844715849, + 0.50101338937930910152, + 0.57558449063515165995, + 0.64317675898520470128, + 0.70355000514714201566, + 0.75669390863372994941, + 0.80279874134324126576, + 0.84221924635075686382, + 0.87543539763040867837, + 0.90301328151357387064, + 0.92556863406861266645, + 0.94373478605275715685, + 0.95813602271021369012, + 0.96936673289691733517, + 0.97797623518666497298, + 0.98445883116743083087, + 0.98924843109013389601, + 0.99271699719682728538, + 0.99517602615532735426, + 0.99688031812819187372, + 0.99803333631543375402, + 0.99879353429880589929, + 0.99928111192179195541, + 0.99958475035151758732, + 0.99976797159956083506, + 0.99987486504878034648, + 0.99993501992508242369, + 0.99996759306794345976, + 0.99998451990227082442, + 0.99999293787666288565, + 0.99999693244919035751, + 0.99999873547186590954, + 0.99999950700571943689, + 0.99999981889371276701, + 0.99999993755407837378, + 0.99999997987450320175, + 0.99999999396413420165, + 0.99999999832336194826, + 0.99999999957078777261, + 0.99999999989927772326, + 0.99999999997845533741, + 0.99999999999582460688, + 0.99999999999927152627, + 0.99999999999988636130 + }, + new[] + { + 0.024539763574649160379, + 0.073525122985671294475, + 0.12222912220155764235, + 0.17046797238201051811, + 0.21806347346971200463, + 0.26484507658344795046, + 0.31065178055284596083, + 0.35533382516507453330, + 0.39875415046723775644, + 0.44078959903390086627, + 0.48133184611690504422, + 0.52028805069123015958, + 0.55758122826077823080, + 0.59315035359195315880, + 0.62695020805104287950, + 0.65895099174335012438, + 0.68913772506166767176, + 0.71750946748732412721, + 0.74407838354734739913, + 0.76886868676824658459, + 0.79191549237614211447, + 0.81326360850297385168, + 0.83296629391941087564, + 0.85108400798784873261, + 0.86768317577564598669, + 0.88283498824466895513, + 0.89661425428007602579, + 0.90909831816302043511, + 0.92036605303195280235, + 0.93049693799715340631, + 0.93957022393327475539, + 0.94766419061515309734, + 0.95485549580502268541, + 0.96121861515111640753, + 0.96682537031235585284, + 0.97174454156548730892, + 0.97604156025657673933, + 0.97977827580061576265, + 0.98301279148110110558, + 0.98579936302528343597, + 0.98818835380074264243, + 0.99022624046752774694, + 0.99195566300267761562, + 0.99341551316926403900, + 0.99464105571251119672, + 0.99566407681695316965, + 0.99651305464025377317, + 0.99721334704346870224, + 0.99778739195890653083, + 0.99825491617199629344, + 0.99863314864067747762, + 0.99893703483351217373, + 0.99917944893488591716, + 0.99937140114093768690, + 0.99952223765121720422, + 0.99963983134560036519, + 0.99973076151980848263, + 0.99980048143113838630, + 0.99985347277311141171, + 0.99989338654759256426, + 0.99992317012928932869, + 0.99994518061445869309, + 0.99996128480785666613, + 0.99997294642523223656, + 0.99998130127012072679, + 0.99998722128200062811, + 0.99999136844834487344, + 0.99999423962761663478, + 0.99999620334716617675, + 0.99999752962380516793, + 0.99999841381096473542, + 0.99999899541068996962, + 0.99999937270733536947, + 0.99999961398855024275, + 0.99999976602333243312, + 0.99999986037121459941, + 0.99999991800479471056, + 0.99999995264266446185, + 0.99999997311323594362, + 0.99999998500307631173, + 0.99999999178645609907, + 0.99999999558563361584, + 0.99999999767323673790, + 0.99999999879798350040, + 0.99999999939177687583, + 0.99999999969875436925, + 0.99999999985405611550, + 0.99999999993088839501, + 0.99999999996803321674, + 0.99999999998556879008, + 0.99999999999364632387, + 0.99999999999727404948, + 0.99999999999886126543, + 0.99999999999953722654, + 0.99999999999981720098, + 0.99999999999992987953 + } + }; + + /// + /// Precomputed weight vector per level. + /// + static readonly double[][] PrecomputedWeights = + new[] + { + new[] + { + 1.5707963267948966192, + 0.230022394514788685, + 0.00026620051375271690866, + 1.3581784274539090834e-12 + }, + new[] + { + 0.96597657941230114801, + 0.018343166989927842087, + 2.1431204556943039358e-7 + }, + new[] + { + 1.3896147592472563229, + 0.53107827542805397476, + 0.076385743570832304188, + 0.0029025177479013135936, + 0.000011983701363170720047, + 1.1631165814255782766e-9 + }, + new[] + { + 1.5232837186347052132, + 1.1934630258491569639, + 0.73743784836154784136, + 0.36046141846934367417, + 0.13742210773316772341, + 0.039175005493600779072, + 0.0077426010260642407123, + 0.00094994680428346871691, + 0.000062482559240744082891, + 1.8263320593710659699e-6, + 1.8687282268736410132e-8, + 4.9378538776631926964e-11 + }, + new[] + { + 1.5587733555333301451, + 1.466014426716965781, + 1.297475750424977998, + 1.0816349854900704074, + 0.85017285645662006895, + 0.63040513516474369106, + 0.44083323627385823707, + 0.290240679312454185, + 0.17932441211072829296, + 0.10343215422333290062, + 0.055289683742240583845, + 0.027133510013712003219, + 0.012083543599157953493, + 0.0048162981439284630173, + 0.0016908739981426396472, + 0.00051339382406790336017, + 0.00013205234125609974879, + 0.000028110164327940134749, + 4.8237182032615502124e-6, + 6.4777566035929719908e-7, + 6.5835185127183396672e-8, + 4.8760060974240625869e-9, + 2.5216347918530148572e-10, + 8.6759314149796046502e-12 + }, + new[] + { + 1.5677814313072218572, + 1.5438811161769592204, + 1.4972262225410362896, + 1.4300083548722996676, + 1.3452788847662516615, + 1.2467012074518577048, + 1.1382722433763053734, + 1.0240449331118114483, + 0.90787937915489531693, + 0.79324270082051671787, + 0.68306851634426375464, + 0.57967810308778764708, + 0.48475809121475539287, + 0.39938474152571713515, + 0.32408253961152890402, + 0.258904639514053516, + 0.20352399885860174519, + 0.15732620348436615027, + 0.11949741128869592428, + 0.089103139240941462841, + 0.065155533432536205042, + 0.046668208054846613644, + 0.032698732726609031113, + 0.022379471063648476483, + 0.014937835096050129696, + 0.0097072237393916892692, + 0.0061300376320830301252, + 0.0037542509774318343023, + 0.0022250827064786427022, + 0.0012733279447082382027, + 0.0007018595156842422708, + 0.00037166693621677760301, + 0.00018856442976700318572, + 0.000091390817490710122732, + 0.000042183183841757600604, + 0.000018481813599879217116, + 7.6595758525203162562e-6, + 2.9916615878138787094e-6, + 1.0968835125901264732e-6, + 3.7595411862360630091e-7, + 1.1992442782902770219e-7, + 3.5434777171421953043e-8, + 9.6498888961089633609e-9, + 2.4091773256475940779e-9, + 5.482835779709497755e-10, + 1.1306055347494680536e-10, + 2.0989335404511469109e-11, + 3.4841937670261059685e-12 + }, + new[] + { + 1.5700420292795931467, + 1.5640214037732320999, + 1.5520531698454121192, + 1.5342817381543034316, + 1.5109197230741697127, + 1.48224329788553807, + 1.4485862549613225916, + 1.4103329714462590129, + 1.3679105116808964881, + 1.3217801174437728579, + 1.2724283455378627082, + 1.2203581095793582207, + 1.1660798699324345766, + 1.1101031939653403796, + 1.0529288799552666556, + 0.99504180404613271514, + 0.93690461274566793366, + 0.87895234555278212039, + 0.82158803526696470334, + 0.7651792989089561367, + 0.71005590120546898385, + 0.65650824613162753076, + 0.60478673057840362158, + 0.55510187800363350959, + 0.5076251588319080997, + 0.4624903980553677613, + 0.41979566844501548066, + 0.37960556938665160999, + 0.3419537959230168323, + 0.30684590941791694932, + 0.27426222968906810637, + 0.24416077786983990868, + 0.21648020911729617038, + 0.19114268413342749532, + 0.16805663794826916233, + 0.14711941325785693248, + 0.12821973363120098675, + 0.11123999898874453035, + 0.096058391865189467849, + 0.082550788110701737654, + 0.070592469906866999352, + 0.060059642358636300319, + 0.05083075757257047107, + 0.042787652157725676034, + 0.035816505604196436523, + 0.029808628117310126969, + 0.024661087314753282511, + 0.020277183817500123926, + 0.016566786254247575375, + 0.013446536605285730674, + 0.010839937168255907211, + 0.0086773307495391815854, + 0.0068957859690660035329, + 0.0054388997976239984331, + 0.0042565295990178580165, + 0.0033044669940348302363, + 0.0025440657675291729678, + 0.0019418357759843675814, + 0.0014690143599429791058, + 0.0011011261134519383862, + 0.00081754101332469493115, + 0.00060103987991147422573, + 0.00043739495615911687786, + 0.00031497209186021200274, + 0.00022435965205008549104, + 0.00015802788400701191949, + 0.00011002112846666697224, + 0.000075683996586201477788, + 0.000051421497447658802092, + 0.0000344921247593431977, + 0.000022832118109036146591, + 0.000014908514031870608449, + 9.5981941283784710776e-6, + 6.0899100320949039256e-6, + 3.8061983264644899045e-6, + 2.3421667208528096843e-6, + 1.4183067155493917523e-6, + 8.4473756384859863469e-7, + 4.9458288702754198508e-7, + 2.8449923659159806339e-7, + 1.6069394579076224911e-7, + 8.9071395140242387124e-8, + 4.8420950198072369669e-8, + 2.579956822953589238e-8, + 1.3464645522302038796e-8, + 6.8784610955899001111e-9, + 3.4371856744650090511e-9, + 1.6788897682161906807e-9, + 8.0099784479729665356e-10, + 3.7299501843052790038e-10, + 1.6939457789411646876e-10, + 7.4967397573818224522e-11, + 3.230446433325236576e-11, + 1.3542512912336274432e-11, + 5.5182369468174885821e-12, + 2.1835922099233609052e-12 + } + }; + + #endregion + + /// + /// Maximum number of iterations, until the asked + /// maximum error is (likely to be) satisfied. + /// + const int NumberOfMaximumLevels = 10; + + /// + /// Abscissa vector per level provider. + /// + IEnumerable _levelAbcissas; + + /// + /// Weight vector per level provider. + /// + IEnumerable _levelWeights; + + /// + /// Approximate the integral by the double exponential transformation + /// + /// The analytic smooth function to integrate. + /// Where the interval starts, inclusive and finite. + /// Where the interval stops, inclusive and finite. + /// The expected relative accuracy of the approximation. + /// Approximation of the finite integral in the given interval. + public double Integrate( + Func f, + double intervalBegin, + double intervalEnd, + double targetRelativeError) + { + if (_levelAbcissas == null) + { + _levelAbcissas = ProvideLevelAbcissas(); + _levelWeights = ProvideLevelWeights(); + } + + return NewtonCotesTrapeziumRule.IntegrateAdaptiveTransformedOdd( + f, + intervalBegin, + intervalEnd, + _levelAbcissas, + _levelWeights, + 1, + targetRelativeError); + } + + /// + /// Abscissa vector per level provider. + /// + /// Level Enumerator. + internal static IEnumerable ProvideLevelAbcissas() + { + for (int i = 0; i < NumberOfMaximumLevels; i++) + { + yield return EvaluateAbcissas(i); + } + } + + /// + /// Weight vector per level provider. + /// + /// Level Enumerator. + internal static IEnumerable ProvideLevelWeights() + { + for (int i = 0; i < NumberOfMaximumLevels; i++) + { + yield return EvaluateWeights(i); + } + } + + /// + /// Compute the abscissa vector for a single level. + /// + /// The level to evaluate the abscissa vector for. + /// Abscissa Vector. + static double[] EvaluateAbcissas(int level) + { + if (level < PrecomputedAbscissas.Length) + { + return PrecomputedAbscissas[level]; + } + + double step = level <= 1 ? 1.0 : (1.0/(2 << (level - 2))); + double offset = level == 0 ? 0.0 : (1.0/(2 << (level - 1))); + int length = level == 0 ? 4 : (3 << (level - 1)); + + double t = 0; + double[] abcissas = new double[length]; + for (int i = 0; i < abcissas.Length; i++) + { + double arg = offset + t; + t += step; + + abcissas[i] = Math.Tanh(Constants.PiOver2*Math.Sinh(arg)); + } + + return abcissas; + } + + /// + /// Compute the weight vector for a single level. + /// + /// The level to evaluate the weight vector for. + /// Weight Vector. + static double[] EvaluateWeights(int level) + { + if (level < PrecomputedWeights.Length) + { + return PrecomputedWeights[level]; + } + + double step = level <= 1 ? 1.0 : (1.0/(2 << (level - 2))); + double offset = level == 0 ? 0.0 : (1.0/(2 << (level - 1))); + int length = level == 0 ? 4 : (3 << (level - 1)); + + double t = 0; + double[] weights = new double[length]; + for (int i = 0; i < weights.Length; i++) + { + double arg = offset + t; + t += step; + + // TODO: reuse abcissas as computed in EvaluateAbcissas + double abcissa = Math.Tanh(Constants.PiOver2*Math.Sinh(arg)); + weights[i] = Constants.PiOver2*(1 - (abcissa*abcissa))*Math.Cosh(arg); + } + + return weights; + } + } +} diff --git a/src/Numerics/Integration/Algorithms/NewtonCotesTrapeziumRule.cs b/src/Numerics/Integration/NewtonCotesTrapeziumRule.cs similarity index 87% rename from src/Numerics/Integration/Algorithms/NewtonCotesTrapeziumRule.cs rename to src/Numerics/Integration/NewtonCotesTrapeziumRule.cs index 6b635b3d..35e4c6fd 100644 --- a/src/Numerics/Integration/Algorithms/NewtonCotesTrapeziumRule.cs +++ b/src/Numerics/Integration/NewtonCotesTrapeziumRule.cs @@ -28,12 +28,12 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.Integration.Algorithms -{ - using System; - using System.Collections.Generic; - using Properties; +using System; +using System.Collections.Generic; +using MathNet.Numerics.Properties; +namespace MathNet.Numerics.Integration +{ /// /// Approximation algorithm for definite integrals by the Trapezium rule of the Newton-Cotes family. /// @@ -59,7 +59,7 @@ namespace MathNet.Numerics.Integration.Algorithms throw new ArgumentNullException("f"); } - return (intervalEnd - intervalBegin) / 2 * (f(intervalBegin) + f(intervalEnd)); + return (intervalEnd - intervalBegin)/2*(f(intervalBegin) + f(intervalEnd)); } /// @@ -86,10 +86,10 @@ namespace MathNet.Numerics.Integration.Algorithms throw new ArgumentOutOfRangeException("numberOfPartitions", Resources.ArgumentPositive); } - double step = (intervalEnd - intervalBegin) / numberOfPartitions; + double step = (intervalEnd - intervalBegin)/numberOfPartitions; double offset = step; - double sum = 0.5 * (f(intervalBegin) + f(intervalEnd)); + double sum = 0.5*(f(intervalBegin) + f(intervalEnd)); for (int i = 0; i < numberOfPartitions - 1; i++) { // NOTE (ruegg, 2009-01-07): Do not combine intervalBegin and offset (numerical stability!) @@ -97,7 +97,7 @@ namespace MathNet.Numerics.Integration.Algorithms offset += step; } - return step * sum; + return step*sum; } /// @@ -121,17 +121,17 @@ namespace MathNet.Numerics.Integration.Algorithms int numberOfPartitions = 1; double step = intervalEnd - intervalBegin; - double sum = 0.5 * step * (f(intervalBegin) + f(intervalEnd)); + double sum = 0.5*step*(f(intervalBegin) + f(intervalEnd)); for (int k = 0; k < 20; k++) { double midpointsum = 0; for (int i = 0; i < numberOfPartitions; i++) { - midpointsum += f(intervalBegin + ((i + 0.5) * step)); + midpointsum += f(intervalBegin + ((i + 0.5)*step)); } midpointsum *= step; - sum = 0.5 * (sum + midpointsum); + sum = 0.5*(sum + midpointsum); step *= 0.5; numberOfPartitions *= 2; @@ -179,9 +179,9 @@ namespace MathNet.Numerics.Integration.Algorithms throw new ArgumentNullException("levelWeights"); } - double linearSlope = 0.5 * (intervalEnd - intervalBegin); - double linearOffset = 0.5 * (intervalEnd + intervalBegin); - targetRelativeError /= 5 * linearSlope; + double linearSlope = 0.5*(intervalEnd - intervalBegin); + double linearOffset = 0.5*(intervalEnd + intervalBegin); + targetRelativeError /= 5*linearSlope; using (var abcissasIterator = levelAbscissas.GetEnumerator()) using (var weightsIterator = levelWeights.GetEnumerator()) @@ -194,10 +194,10 @@ namespace MathNet.Numerics.Integration.Algorithms double[] abcissasL1 = abcissasIterator.Current; double[] weightsL1 = weightsIterator.Current; - double sum = f(linearOffset) * weightsL1[0]; + double sum = f(linearOffset)*weightsL1[0]; for (int i = 1; i < abcissasL1.Length; i++) { - sum += weightsL1[i] * (f((linearSlope * abcissasL1[i]) + linearOffset) + f(-(linearSlope * abcissasL1[i]) + linearOffset)); + sum += weightsL1[i]*(f((linearSlope*abcissasL1[i]) + linearOffset) + f(-(linearSlope*abcissasL1[i]) + linearOffset)); } sum *= step; @@ -212,11 +212,11 @@ namespace MathNet.Numerics.Integration.Algorithms double midpointsum = 0; for (int i = 0; i < abcissas.Length; i++) { - midpointsum += weights[i] * (f((linearSlope * abcissas[i]) + linearOffset) + f(-(linearSlope * abcissas[i]) + linearOffset)); + midpointsum += weights[i]*(f((linearSlope*abcissas[i]) + linearOffset) + f(-(linearSlope*abcissas[i]) + linearOffset)); } midpointsum *= step; - sum = 0.5 * (sum + midpointsum); + sum = 0.5*(sum + midpointsum); step *= 0.5; double delta = Math.Abs(sum - midpointsum); @@ -227,7 +227,7 @@ namespace MathNet.Numerics.Integration.Algorithms continue; } - double r = Math.Log(delta) / Math.Log(previousDelta); + double r = Math.Log(delta)/Math.Log(previousDelta); previousDelta = delta; if (r > 1.9 && r < 2.1) @@ -242,7 +242,7 @@ namespace MathNet.Numerics.Integration.Algorithms } } - return sum * linearSlope; + return sum*linearSlope; } } } diff --git a/src/Numerics/Integration/Algorithms/SimpsonRule.cs b/src/Numerics/Integration/SimpsonRule.cs similarity index 88% rename from src/Numerics/Integration/Algorithms/SimpsonRule.cs rename to src/Numerics/Integration/SimpsonRule.cs index d4fcf35d..721cf37f 100644 --- a/src/Numerics/Integration/Algorithms/SimpsonRule.cs +++ b/src/Numerics/Integration/SimpsonRule.cs @@ -28,12 +28,12 @@ // OTHER DEALINGS IN THE SOFTWARE. // -namespace MathNet.Numerics.Integration.Algorithms -{ - using System; - using NumberTheory; - using Properties; +using System; +using MathNet.Numerics.NumberTheory; +using MathNet.Numerics.Properties; +namespace MathNet.Numerics.Integration +{ /// /// Approximation algorithm for definite integrals by Simpson's rule. /// @@ -56,8 +56,8 @@ namespace MathNet.Numerics.Integration.Algorithms throw new ArgumentNullException("f"); } - double midpoint = (intervalEnd + intervalBegin) / 2; - return (intervalEnd - intervalBegin) / 6 * (f(intervalBegin) + f(intervalEnd) + (4 * f(midpoint))); + double midpoint = (intervalEnd + intervalBegin)/2; + return (intervalEnd - intervalBegin)/6*(f(intervalBegin) + f(intervalEnd) + (4*f(midpoint))); } /// @@ -89,8 +89,8 @@ namespace MathNet.Numerics.Integration.Algorithms throw new ArgumentException(Resources.ArgumentEven, "numberOfPartitions"); } - double step = (intervalEnd - intervalBegin) / numberOfPartitions; - double factor = step / 3; + double step = (intervalEnd - intervalBegin)/numberOfPartitions; + double factor = step/3; double offset = step; int m = 4; @@ -98,12 +98,12 @@ namespace MathNet.Numerics.Integration.Algorithms for (int i = 0; i < numberOfPartitions - 1; i++) { // NOTE (cdrnet, 2009-01-07): Do not combine intervalBegin and offset (numerical stability) - sum += m * f(intervalBegin + offset); + sum += m*f(intervalBegin + offset); m = 6 - m; offset += step; } - return factor * sum; + return factor*sum; } } } diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index a333bc40..39987e82 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -392,10 +392,10 @@ - - - - + + + + diff --git a/src/UnitTests/IntegrationTests/IntegrationTest.cs b/src/UnitTests/IntegrationTests/IntegrationTest.cs index e006c095..d2ed8083 100644 --- a/src/UnitTests/IntegrationTests/IntegrationTest.cs +++ b/src/UnitTests/IntegrationTests/IntegrationTest.cs @@ -28,7 +28,6 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests { using System; using Integration; - using Integration.Algorithms; using NUnit.Framework; ///