C# :: Aufgabe #195

3 Lösungen Lösungen öffentlich

Berechnung des mittleren Punktabstandes in geometrischen Formen

Anfänger - C# von hollst - 17.12.2017 um 21:21 Uhr
In mathebord regt man sich u. a. darüber auf,
dass man bei einer Monte-Carlo-Simulation zur Berechnung des mittleren Abstandes
zweier Punkte in einem Rechteck über 4 Sekunden benötigt (bei 1 Mio. zufälligen Experimenten).
Außerdem wird in der Antwort vom 22.10.2013 (23:10) auch noch ein falsches Ergebnis unkommentiert präsentiert.

Nun gut, das war vor etwa 5 Jahren und Hard- sowie verwendete Software waren auch nicht Stand der damaligen Technik.
Wir wollen versuchen, es hier schneller und besser zu bewerkstelligen.

Aufgabenstellung: Für die fünf (euklidischen) Formen gerade Linie, Quadrat, Rechteck, Kreis und Ellipse
schätze man den mittleren Abstand zweier Punkte innerhalb der genannten geometrischen Formen ab, bezogen auf den Umfang der entsprechenden Form. Die Umfangspunkte gehören mit zu der Form. Für die Linie sei der Umfang die doppelte Linienlänge
(Rechteck mit zwei verschwindend kleinen Parallelseiten).

Die Simulation soll mittels 1.000.000 Zufallsexperimenten vorgenommen werden.

Lösungen:

vote_ok
von hollst (12440 Punkte) - 22.12.2017 um 11:33 Uhr
Quellcode ausblenden C#-Code
using System;
using static System.Console;
using static System.ConsoleKey;

//http://www.matheboard.de/archive/529571/thread.html
//http://www.matheboard.de/archive/25782/thread.html

namespace mittlerer_abstand_zweier_kreispunkte
{
    static class Program
    {
        static string NL = Environment.NewLine, Version = "version 18.12.2017" + NL;

        static void Main()
        {
            Version.MessageLine();
            bool bo_loop = true;
            while (bo_loop)
            {
                int max_trails = (int)1e+6;
                double r1 = 0, r2 = 0;

                bool bo_input_okay = false;
                while (!bo_input_okay)
                {
                    "give me r1: ".Message();
                    string string_input = ReadLine();
                    bo_input_okay = double.TryParse(string_input, out r1) && (r1 >= 0.0);
                }
                bo_input_okay = false;
                while (!bo_input_okay)
                {
                    "give me r2: ".Message();
                    string string_input = ReadLine();
                    bo_input_okay = double.TryParse(string_input, out r2) && (r2 >= 0.0);
                }

                string[] info = { "line     ", "quad     ", "rect     ", "circle   ", "ellipse  " };

                double[][] radii = new double[][]  //  die formenmittelpunkte liegen bei (0.0, 0.0)
                {
                    new double[]{r1, 0.0},                    new double[]{r2, 0.0},    //line

                    new double[]{r1,  r1},                    new double[]{r2,  r2},    //quadrat
                    new double[]{r1,  r2},                    new double[]{r2,  r1},    //recheck

                    new double[]{r1,  r1},                    new double[]{r2,  r2},    //kreis
                    new double[]{r1,  r2},                    new double[]{r2,  r1}     //ellipse
                };

                Delegate [][] deligierte = new Delegate[][] //  regarding radii
                {
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Quad_Rect_Line, (IsIn_Area)MyDistance_Delegates.IsIn_Quad_Rect_Line },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Quad_Rect_Line, (IsIn_Area)MyDistance_Delegates.IsIn_Quad_Rect_Line },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Quad_Rect_Line, (IsIn_Area)MyDistance_Delegates.IsIn_Quad_Rect_Line },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Quad_Rect_Line, (IsIn_Area)MyDistance_Delegates.IsIn_Quad_Rect_Line },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Quad_Rect_Line, (IsIn_Area)MyDistance_Delegates.IsIn_Quad_Rect_Line },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Quad_Rect_Line, (IsIn_Area)MyDistance_Delegates.IsIn_Quad_Rect_Line },

                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Circle_Ellipse, (IsIn_Area)MyDistance_Delegates.IsIn_Circle_Ellipse },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Circle_Ellipse, (IsIn_Area)MyDistance_Delegates.IsIn_Circle_Ellipse },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Circle_Ellipse, (IsIn_Area)MyDistance_Delegates.IsIn_Circle_Ellipse },
                    new Delegate[] { (Circumference)MyDistance_Delegates.Circumference_Circle_Ellipse, (IsIn_Area)MyDistance_Delegates.IsIn_Circle_Ellipse }
                };

                (NL + "Results from Class_Mean_Distances in a Loop").MessageLine();
                for (var i = 0; i < radii.Length; i++)
                {
                    Class_Mean_Distances
                        cmd = new Class_Mean_Distances(radii[i][0], radii[i][1], max_trails, (Circumference)deligierte[i][0], (IsIn_Area)deligierte[i][1]);
                    $"{info[i / 2]}max_trails: {max_trails.ToString("n0")}    r1: {radii[i][0],10}   result = {(100.0 * cmd.result).ToString("0.0000"),5} %    r2: {radii[i][1],10}".MessageLine();
                }

                (NL + "test circumference calculation for ellipse").MessageLine();
                Class_Mean_Distances
                    cmd_exact = new Class_Mean_Distances(r1, r2, max_trails, MyDistance_Delegates.Circumference_Ellipse_exact, MyDistance_Delegates.IsIn_Circle_Ellipse);
                $"exact:   max_trails: {max_trails.ToString("n0")}    r1: {r1,10}   result = {(100.0 * cmd_exact.result).ToString("0.0000"),5} %    r2: {r2,10}   U_exact: {cmd_exact.U.ToString()}".MessageLine();
                    cmd_exact = new Class_Mean_Distances(r1, r2, max_trails, MyDistance_Delegates.Circumference_Circle_Ellipse, MyDistance_Delegates.IsIn_Circle_Ellipse);
                $"approx:  max_trails: {max_trails.ToString("n0")}    r1: {r1,10}   result = {(100.0 * cmd_exact.result).ToString("0.0000"),5} %    r2: {r2,10}   U      : {cmd_exact.U.ToString()}".MessageLine();

                ConsoleKeyInfo ki = "press any key (exit => ESC)".ReadKey(bo_NL: true);
                bo_loop = !(ki.Key == Escape);
            }
            "ready, press a key".ReadKey(true);
        }

        #region expansions

        public static void Message(this string s) => Write(s);

        public static void MessageLine(this string s) => WriteLine(s);

        public static ConsoleKeyInfo ReadKey(this string s, bool bo_NL = false)
        {
            if (bo_NL) s.MessageLine(); else s.Message();
            return System.Console.ReadKey(true);
        }
        #endregion
    }

    // define delegates

    public delegate double Circumference(double r1, double r2);
    public delegate bool IsIn_Area(double r1, double r2, double x, double y);    

    public static class MyDistance_Delegates // create methods for the delegates
    {
        public static bool IsIn_Quad_Rect_Line(double r1, double r2, double x, double y) => true;   //quad: r2 = r1; line: r2 = 0

        public static double Circumference_Quad_Rect_Line(double r1, double r2) => 4 * (r1 + r2);   //quad: r2 = r1; line: r2 = 0

        public static bool IsIn_Circle_Ellipse(double r1, double r2, double x, double y)            //circle: r2 = r1
            => ((x - r1) * (x - r1) * r2 * r2 + (y - r2) * (y - r2) * r1 * r1 <= r1 * r1 * r2 * r2);

        public static double Circumference_Circle_Ellipse(double r1, double r2)                     //circle: r2 = r1
        {
            double lambda = (r1 - r2) / (r1 + r2);  //approximation from Ramanujan; for circles (lambda = 0) exact
            return Math.PI * (r1 + r2) * (1.0 + 3 * lambda * lambda / (10.0 + Math.Sqrt(4.0 - 3.0 * lambda * lambda)));
        }

        public static double Circumference_Ellipse_exact(double r1, double r2)                     //circle: r2 = r1
        //see https://de.wikipedia.org/wiki/Ellipse
        {
            double lambda = (r1 - r2) / (r1 + r2), lambda_quadrat = lambda * lambda;
            double error = 1.0E-16; //1.0E-38;

            double sum = 1.0, sum_old = 0.0, sum_n = 1.0;
            double no = 1, nu = 1;
            int n = 0;
            while((sum - sum_old) > error)
            {
                sum_old = sum;
                no = (2 * n - 1);
                nu = (2 * n + 2);
                double q = (double)no / (double)nu;
                sum_n *= lambda_quadrat * q * q;
                sum += sum_n;
                n++;
            }
$"n from Circumference_Ellipse_exact: {n.ToString("n0")}".MessageLine();
            return sum * Math.PI * (r1 + r2);
        }
    }

    public class Class_Mean_Distances   //  with use of delegates
    {
        private double r1, r2;
        private int max_loops;

        private Circumference circumstance;
        private IsIn_Area bo_isinarea;

        private Random rand = new Random();

        public double U, result;

        public Class_Mean_Distances(double r1, double r2, int max_loops, Circumference circumstance, IsIn_Area bo_isinarea)
        {
            this.r1 = r1;
            this.r2 = r2;
            this.max_loops = max_loops;
            this.circumstance = circumstance;
            this.bo_isinarea = bo_isinarea;

            this.run();
        }

        private void run()
        {
            int counter = 0;
            double distance = 0.0;
            while (counter < max_loops)
            {
                double[][] xy = new double[2][];
                for(var i = 0; i < xy.Length; i++)
                {
                    bool bo = false;
                    while (!bo)
                    {
                        double x = rand.NextDouble() * 2 * r1;
                        double y = rand.NextDouble() * 2 * r2;
                        xy[i] = new double[] { x, y };
                        bo = bo_isinarea(r1, r2, xy[i][0], xy[i][1]);
                    }
                }
                double act_distance = Math.Sqrt((xy[0][0] - xy[1][0]) * (xy[0][0] - xy[1][0]) + 
                                                (xy[0][1] - xy[1][1]) * (xy[0][1] - xy[1][1]));
                distance += act_distance;
                counter++;
            }
            U = circumstance(r1, r2);
            this.result = (distance / max_loops) / U;
        }
    }
}
vote_ok
von daniel59 (4260 Punkte) - 11.01.2018 um 09:46 Uhr
Quellcode ausblenden C#-Code
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Threading.Tasks;

namespace ConsoleMittlererPunktabstand
{
    class Program
    {
        static readonly Random rnd = new Random();
        static void Main(string[] args)
        {
            int iterations = 1000000;
            double maxWidth = 1;
            double maxHeight = 1;

            Stopwatch sw = new Stopwatch();
            sw.Start();
            Task<double> taskLineExperiment = Task.Factory.StartNew(() => LineExperiment(iterations, maxWidth));
            Task<double> taskSquareExperiment = Task.Factory.StartNew(() => SquareExperiment(iterations, maxWidth));
            Task<double> taskRectangleExperiment = Task.Factory.StartNew(() => RectangleExperiment(iterations, maxWidth, maxHeight));
            Task<double> taskCircleExperiment = Task.Factory.StartNew(() => CircleExperiment(iterations, maxWidth));
            Task<double> taskEllipseExperiment = Task.Factory.StartNew(() => EllipseExperiment(iterations, maxWidth, maxHeight));

            Task.WaitAll(taskLineExperiment, taskSquareExperiment, taskRectangleExperiment, taskCircleExperiment, taskEllipseExperiment);
            sw.Stop();


            Console.WriteLine($"{nameof(LineExperiment)}: {taskLineExperiment.Result}");
            Console.WriteLine($"{nameof(SquareExperiment)}: {taskSquareExperiment.Result}");
            Console.WriteLine($"{nameof(RectangleExperiment)}: {taskRectangleExperiment.Result}");
            Console.WriteLine($"{nameof(CircleExperiment)}: {taskCircleExperiment.Result}");
            Console.WriteLine($"{nameof(EllipseExperiment)}: {taskEllipseExperiment.Result}");
            Console.WriteLine($"Benötigte Zeit: {sw.Elapsed}");
            Console.ReadLine();
        }

        static double LineExperiment(int iterations, double maxWidth)
        {
            return CreateManyForms(EuclidianType.Line, iterations, maxWidth).Average(a => CalcDistance(a));
        }

        static double SquareExperiment(int iterations, double maxWidth)
        {
            return CreateManyForms(EuclidianType.Square, iterations, maxWidth).Average(a => CalcDistance(a));
        }

        static double RectangleExperiment(int iterations, double maxWidth, double maxHeight)
        {
            return CreateManyForms(EuclidianType.Rectangle, iterations, maxWidth, maxHeight).Average(a => CalcDistance(a));
        }

        static double CircleExperiment(int iterations, double maxWidth)
        {
            return CreateManyForms(EuclidianType.Circle, iterations, maxWidth).Average(a => CalcDistance(a));
        }

        static double EllipseExperiment(int iterations, double maxWidth, double maxHeight)
        {
            return CreateManyForms(EuclidianType.Ellipse, iterations, maxWidth, maxHeight).Average(a => CalcDistance(a));
        }

        static double CalcDistance(EuclidianForm form)
        {
            if (form != null)
            {
                Point first = form.GeneratePoint();
                Point second = form.GeneratePoint();

                Vector diff = second - first;

                return diff.Length / form.Umfang;
            }
            else
            {
                return double.NaN;
            }
        }

        static EuclidianForm CreateRandomForm(EuclidianType type, double maxWidth, double maxHeight = 0)
        {
            EuclidianForm form = new EuclidianForm();
            form.Type = type;
            lock (rnd)
            {
                double w = rnd.NextDouble();
                double h;
                form.Width = maxWidth * w;
                if (maxHeight != 0)
                {
                    h = rnd.NextDouble();
                    form.Height = maxHeight * h;
                }
            }
            return form;
        }

        static IEnumerable<EuclidianForm> CreateManyForms(EuclidianType type, int iterations, double maxWidth, double maxHeight = 0)
        {
            for (int i = 0; i < iterations; i++)
            {
                yield return CreateRandomForm(type, maxWidth, maxHeight);
            }
        }
    }

    class EuclidianForm
    {
        private static readonly Random rnd = new Random();

        public EuclidianType Type { get; set; }

        private double width = 1;
        public double Width
        {
            get { return width; }
            set
            {
                width = value;
                if (Type == EuclidianType.Circle || Type == EuclidianType.Square)
                { height = value; }
            }
        }

        private double height = 1;
        public double Height
        {
            get { return height; }
            set
            {
                if (Type != EuclidianType.Line)
                {
                    height = value;
                    if (Type == EuclidianType.Circle || Type == EuclidianType.Square)
                    { width = value; }
                }
            }
        }

        public double Umfang => calcUmfang();


        #region Erzeugen von zufälligen Punkten
        public Point GeneratePoint()
        {
            switch (Type)
            {
                case EuclidianType.Line:
                    return GeneratePointLine();

                case EuclidianType.Square:
                    return GeneratePointSquare();

                case EuclidianType.Rectangle:
                    return GeneratePointRectangle();

                case EuclidianType.Circle:
                    return GeneratePointCircle();

                case EuclidianType.Ellipse:
                    return GeneratePointEllipse();

                default:
                    throw new ArgumentException();
            }
        }

        private Point GeneratePointLine()
        {
            Point p = new Point();
            lock (rnd)
            {
                p.X = Width * rnd.NextDouble();
                p.Y = 0;
            }
            return p;
        }

        private Point GeneratePointSquare()
        {
            Point p = new Point();
            lock (rnd)
            {
                p.X = Width * rnd.NextDouble();
                p.Y = Width * rnd.NextDouble();
            }
            return p;
        }

        private Point GeneratePointRectangle()
        {
            Point p = new Point();
            lock (rnd)
            {
                p.X = Width * rnd.NextDouble();
                p.Y = Height * rnd.NextDouble();
            }
            return p;
        }

        private Point GeneratePointCircle()
        {
            Point p = new Point();
            lock (rnd)
            {
                double radius = Width / 2;

                double x = rnd.NextDouble() * Width - radius;

                double maxY = Math.Sqrt(radius * radius - x * x);
                double minY = -maxY;

                double y = (maxY - minY) * rnd.NextDouble() + minY;

                p.X = x;
                p.Y = y;
            }

            return p;
        }

        private Point GeneratePointEllipse()
        {
            Point p = new Point();

            lock (rnd)
            {
                double a = Width / 2;
                double b = Height / 2;

                double x = rnd.NextDouble() * Width - a;
                double maxY = Math.Sqrt((b * b) - ((b * b) / (a * a)) * (x * x));

                double minY = -maxY;

                double y = (maxY - minY) * rnd.NextDouble() + minY;

                p.X = x;
                p.Y = y;
            }

            return p;
        }
        #endregion

        private double calcUmfang()
        {
            switch (Type)
            {
                case EuclidianType.Line:
                    return 2 * Width;

                case EuclidianType.Square:
                    return 4 * Width;

                case EuclidianType.Rectangle:
                    return 2 * Width + 2 * Height;

                case EuclidianType.Circle:
                    return Math.PI * Width;

                case EuclidianType.Ellipse:
                    //Näherungsformel nach Ramanujan
                    double a = Width / 2;
                    double b = Height / 2;
                    double l = Math.Pow(((a - b) / (a + b)), 2) * 3;
                    return Math.PI * (a + b) * (1 + (l / (10 + Math.Sqrt(4 - l))));

                default:
                    return double.NaN;
            }
        }


    }

    class Point
    {
        public double X { get; set; }
        public double Y { get; set; }

        public static Vector operator +(Point left, Point right) => new Vector() { X = left.X + right.X, Y = left.Y + right.Y };
        public static Vector operator -(Point left, Point right) => new Vector() { X = left.X - right.X, Y = left.Y - right.Y };
    }

    class Vector : Point
    {
        public double Length => Math.Sqrt(X * X + Y * Y);
    }

    enum EuclidianType
    {
        Line, Square, Rectangle, Circle, Ellipse
    }
}
vote_ok
von luckyman81 (550 Punkte) - 23.04.2020 um 22:54 Uhr
Quellcode ausblenden C#-Code
using System;
using System.Collections.Generic;
using System.Linq;

namespace CS_Aufgabe_195_AverageDistance
{
    class Program
    {
        static void Main(string[] args)
        {
            Random rnd = new Random();
            List<double> dist = new List<double>();
            const double a = 1e0; // Seiten des Rechtecks bzw große und kleine Halbachse
            double b = rnd.NextDouble();
            double iter = 1e6;

            // Straight line
            for (int i = 0; i < iter; i++)
            {
                double x = rnd.NextDouble();
                double y = rnd.NextDouble();
                dist.Add(Math.Abs(x - y));
            }
            Console.WriteLine($"Line:   {dist.Average()}");
            dist.Clear();

            // Square
            for (int i = 0; i < iter; i++)
            {
                double x1 = rnd.NextDouble();
                double x2 = rnd.NextDouble();
                double y1 = rnd.NextDouble();
                double y2 = rnd.NextDouble();
                dist.Add(Math.Sqrt(Math.Pow(x1 - x2, 2e0) + Math.Pow(y1 - y2, 2e0)));
            }
            Console.WriteLine($"Square:   {dist.Average()}");
            dist.Clear();

            // Rectangle
            for (int i = 0; i < iter; i++)
            {
                double x1 = rnd.NextDouble() * a;
                double x2 = rnd.NextDouble() * a;
                double y1 = rnd.NextDouble() * b;
                double y2 = rnd.NextDouble() * b;
                dist.Add(Math.Sqrt(Math.Pow(x1 - x2, 2e0) + Math.Pow(y1 - y2, 2e0)));
            }
            Console.WriteLine($"Rectangle:   {dist.Average()}");
            dist.Clear();

            // Circle
            for (int i = 0; i < iter; i++)
            {
                double ra = rnd.NextDouble();
                double rb = rnd.NextDouble();
                double thetaa = rnd.NextDouble() * Math.PI;
                double thetab = rnd.NextDouble() * Math.PI;
                dist.Add(Math.Sqrt(Math.Pow(ra, 2e0) + Math.Pow(rb, 2e0) - 2e0 * ra * rb * Math.Cos(Math.Abs(thetaa-thetab))));

            }
            Console.WriteLine($"Circle:   {dist.Average()}");
            dist.Clear();

            // Ellipse
            for (int i = 0; i < iter; i++)
            {
                double thetaa = rnd.NextDouble() * Math.PI;
                double rthetaa = a * b / Math.Sqrt(a * a * Math.Pow(Math.Sin(thetaa), 2e0) + b * b * Math.Pow(Math.Cos(thetaa), 2e0));

                double thetab = rnd.NextDouble() * Math.PI;
                double rthetab = a * b / Math.Sqrt(a * a * Math.Pow(Math.Sin(thetab), 2e0) + b * b * Math.Pow(Math.Cos(thetab), 2e0));
                
                dist.Add(Math.Sqrt(Math.Pow(rthetaa, 2e0) + Math.Pow(rthetab, 2e0) - 2e0 * rthetaa * rthetab * Math.Cos(Math.Abs(thetaa-thetab))));

            }
            Console.WriteLine($"Ellipse:   {dist.Average()}");
            dist.Clear();

            Console.ReadKey();
        }
    }
}