Udregne en "improper integral" i C#

1) formlen er forkert i jeres PNG
2) dette er en meget kendt problem stilling med nogle meget avancerede loesninger
3) men p.g.a. tidsmangel faar I kun er meget simpel loesning

using System;

namespace E
{
    public class Calculus
    {
        private const int N = 1000000;
        public delegate double Func(double x);
        public static double Integrate(Func f, double start, double finish)
        {
            double eps = (finish - start) / N;
            double sum = 0;
            for(int i = 0; i < N; i++)
            {
                sum += eps*f(start + eps * (0.5 + i));
            }
            return sum;
        }
    }
    public class Program
    {
        public static double StandardNormal(double x)
        {
            return Math.Exp(-x*x/2)/Math.Sqrt(2*Math.PI);
        }
        public static void Main(string[] args)
        {
            //Console.WriteLine(Math.PI / 2);
            //Console.WriteLine(Calculus.Integrate((double x) => Math.Sqrt(1 - x*x), -1, 1));
            Console.WriteLine(Calculus.Integrate(StandardNormal, -10, 1.345067));
            Console.ReadKey();
        }
    }
}

Hvis jeg faar tid i weekende kan I faa en mere avanceret loesning.

Mange tak arne_v :-) Vi vil prøve ovenstående og vi vil være virkelig taknemmelige hvis du gad bruge tid på en mere avanceret løsning. Specielt da vi gerne vil have det så præcist som muligt og vi føler vi har stødt panden mod en mur, men som du nu har gjort væsenligt mindre.

Den er skam rimeligt præcis. I har næppe data som har brug for større præcision.

Forbedringerne vil primært få den til at køre hurtigere. Ovenstående deler kurven op i en million små rektanker og beregner arealet af det.

Noget mere avanceret kode:

using System;

namespace E
{
    public class Calculus
    {
        public delegate double Func(double x);
        public static double RectangularIntegrate(Func f, double start, double finish, int nintval)
        {
            double eps = (finish - start) / nintval;
            double sum = 0;
            for(int i = 0; i < nintval; i++)
            {
                sum += f(start + eps * (0.5 + i));
            }
            return eps * sum;
        }
        public static double TrapezIntegrate(Func f, double start, double finish, int nintval)
        {
            double eps = (finish - start) / nintval;
            double sum = 0.5 * f(start);
            for(int i = 1; i < nintval; i++)
            {
                sum += f(start + eps * i);
            }
            sum += 0.5 * f(finish);
            return eps * sum;
        }
        public static double SimpsonIntegrate(Func f, double start, double finish, int nintval)
        {
            double eps = (finish - start) / (2 * nintval);
            double sum = f(start);
            for(int i = 1; i <= nintval; i++)
            {
                sum += 4 * f(start + eps * (2 * i - 1));
            }
            for(int i = 1; i < nintval; i++)
            {
                sum += 2 * f(start + eps * 2 * i);
            }
            sum += f(finish);
            return eps * sum / 3;
        }
    }
    public class Program
    {
        private const string FMT = "{0,-30} : {1:f12} {2}";
        public static double StandardNormal(double x)
        {
            return Math.Exp(-x*x/2)/Math.Sqrt(2*Math.PI);
        }
        private static double[] poly1 = {  0.000124818087,
                                          -0.001075204047,
                                          0.005198775019,
                                          -0.019198292004,
                                          0.059054035642,
                                          -0.151968751364,
                                          0.319152932604,
                                          -0.531923007300,
                                          0.797884560593 };
        private static double[] poly2 = { -0.000045255659,
                                          0.000152529290,
                                          -0.000019538132,
                                          -0.000676904986,
                                          0.001390604284,
                                          -0.000794620820,
                                          -0.002034254874,
                                          0.006549791214,
                                          -0.010557625006,
                                          0.011630447319,
                                          -0.009279453341,
                                          0.005353579108,
                                          -0.002141268741,
                                          0.000535310849,
                                          0.999936657524 };
        public static double StandardNormalAccAlgo1(double x)
        {
            double tmp;
            double absx = Math.Abs(x);
            if(absx == 0) {
              tmp = 0;
            } else if(absx > 6) {
              tmp = 1;
            } else if(absx < 2) {
              double v = absx*absx/4;
              tmp = 0;
              for(int i = 0; i < poly1.Length; i++) tmp = tmp * v + poly1[i];
              tmp = tmp * absx;
            } else {
              double v = absx / 2 - 2;
              tmp = 0;
              for(int i = 0; i < poly2.Length; i++) tmp = tmp * v + poly2[i];
            }
            if(x > 0)
              return (tmp+1)/2;
            else
              return (1-tmp)/2;
        }
        private const double b1 =  0.319381530;
        private const double b2 = -0.356563782;
        private const double b3 =  1.781477937;
        private const double b4 = -1.821255978;
        private const double b5 =  1.330274429;
        private const double p  =  0.2316419;
        private const double c  =  0.39894228;
        public static double StandardNormalAccAlgo2(double x)
        {
            if(x >= 0.0)
            {
                double t = 1.0 / ( 1.0 + p * x );
                return (1.0 - c * Math.Exp( -x * x / 2.0 ) * t * ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
            }
            else
            {
                double t = 1.0 / ( 1.0 - p * x );
                return ( c * Math.Exp( -x * x / 2.0 ) * t * ( t *( t * ( t * ( t * b5 + b4 ) + b3 ) + b2 ) + b1 ));
            }
        }
        private static string Elapsed(DateTime dt1, DateTime dt2, int rep)
        {
            return string.Format("({0:F8} seconds)", (dt2-dt1).TotalSeconds/rep);
        }
        private const int REP = 100;
        private static void TestRectangular(int n)
        {
            double res = 0;
            DateTime dt1 = DateTime.Now;
            for(int i = 0; i < REP; i++)
            {
                res = Calculus.RectangularIntegrate(StandardNormal, -1, 1, n);
            }
            DateTime dt2 = DateTime.Now;
            Console.WriteLine(FMT, "Rectangular " + n + " slices", res, Elapsed(dt1, dt2, REP));
        }
        private static void TestTrapez(int n)
        {
            double res = 0;
            DateTime dt1 = DateTime.Now;
            for(int i = 0; i < REP; i++)
            {
                res = Calculus.TrapezIntegrate(StandardNormal, -1, 1, n);
            }
            DateTime dt2 = DateTime.Now;
            Console.WriteLine(FMT, "Trapez " + n + " slices", res, Elapsed(dt1, dt2, REP));
        }
        private static void TestSimpson(int n)
        {
            double res = 0;
            DateTime dt1 = DateTime.Now;
            for(int i = 0; i < REP; i++)
            {
                res = Calculus.SimpsonIntegrate(StandardNormal, -1, 1, n);
            }
            DateTime dt2 = DateTime.Now;
            Console.WriteLine(FMT, "Simpson " + n + " slices", res, Elapsed(dt1, dt2, REP));
        }
        private static void TestAlgo1()
        {
            double res = 0;
            DateTime dt1 = DateTime.Now;
            for(int i = 0; i < REP; i++)
            {
                res = StandardNormalAccAlgo1(1) - StandardNormalAccAlgo1(-1);
            }
            DateTime dt2 = DateTime.Now;
            Console.WriteLine(FMT, "Algorithm 1", res, Elapsed(dt1, dt2, REP));
        }
        private static void TestAlgo2()
        {
            double res = 0;
            DateTime dt1 = DateTime.Now;
            for(int i = 0; i < REP; i++)
            {
                res = StandardNormalAccAlgo2(1) - StandardNormalAccAlgo1(-1);
            }
            DateTime dt2 = DateTime.Now;
            Console.WriteLine(FMT, "Algorithm 2", res, Elapsed(dt1, dt2, REP));
        }
        public static void Main(string[] args)
        {
            Console.WriteLine("Test");
            Console.WriteLine("----");
            Console.WriteLine(FMT, "Correct", 0.682689492137, "");
            TestRectangular(10);
            TestRectangular(100);
            TestRectangular(1000);
            TestRectangular(10000);
            TestRectangular(100000);
            TestRectangular(1000000);
            TestTrapez(10);
            TestTrapez(100);
            TestTrapez(1000);
            TestTrapez(10000);
            TestTrapez(100000);
            TestTrapez(1000000);
            TestSimpson(10);
            TestSimpson(100);
            TestSimpson(1000);
            TestSimpson(10000);
            TestSimpson(100000);
            TestSimpson(1000000);
            TestAlgo1();
            TestAlgo2();
            Console.WriteLine("Result");
            Console.WriteLine("------");
            Console.WriteLine(FMT, "Rectangular 1000 slices", Calculus.RectangularIntegrate(StandardNormal, -10, 1.345067, 1000), "");
            Console.WriteLine(FMT, "Trapez 1000 slices", Calculus.TrapezIntegrate(StandardNormal, -10, 1.345067, 1000), "");
            Console.WriteLine(FMT, "Simpson 1000 slices", Calculus.SimpsonIntegrate(StandardNormal, -10, 1.345067, 1000), "");
            Console.WriteLine(FMT, "Rectangulat 1000000 slices", Calculus.RectangularIntegrate(StandardNormal, -10, 1.345067, 1000000), "");
            Console.WriteLine(FMT, "Trapez 1000000 slices", Calculus.TrapezIntegrate(StandardNormal, -10, 1.345067, 1000000), "");
            Console.WriteLine(FMT, "Simpson 1000000 slices", Calculus.SimpsonIntegrate(StandardNormal, -10, 1.345067, 1000000), "");
            Console.WriteLine(FMT, "Algorithm 1", StandardNormalAccAlgo1(1.345067), "");
            Console.WriteLine(FMT, "Algorithm 2", StandardNormalAccAlgo2(1.345067), "");
            Console.ReadKey();
        }
    }
}

OK ?