Численное Интегрирование

МОСКОВСКИЙ ЭНЕРГЕТИЧЕСКИЙ ИНСТИТУТ

(НАУЧНО-ИССЛЕДОВАТЕЛЬСКИЙ ИНСТИТУТ)

 

 

 

 

 

 

 

 

 

 

 

 

Лабораторная  работа № 8

 «Численное интегрирование».

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Группа:                     А -13 -12

Студент:             Ефанов Андрей

Вариант № 11.

 

 

 

 

Москва 2014

Задача 8.1. . Вычислить значение интеграла , где , с помощью квадратурных формул трапеций и Симпсона с точностью . Предварительно оценить

шаг интегрирования, при котором достигается заданная точность. Сравнить время вычисления интеграла.

 

Теоретический материал.

Наиболее широко для вычисления интегралов вида на практике используют квадратурные формулы – приближенные равенства вида:

 

узлы квадратурной формулы, – веса квадратурной формулы.

 

Можно вывести простейшие квадратурные формулы, исходя из наглядных геометрических соображений. Будем интерпретировать интеграл как площадь криволинейной трапеции.

Разобьем отрезок [a,b] на элементарные отрезки [xi-1,xi]. Интеграл разобьется на сумму элементарных интегралов Ii Обозначим , i = 1..n.

 

Формула трапеций

Соединяем отрезком точки ( и (. Получим формулу трапеций:

 

Остаточный член формулы трапеций:

 

 

Формула Симпсона

Если площадь элементарной криволинейной трапеции заменить площадью фигуры, расположенной под параболой, проходящей через точки , то получим приближенное равенство: , – интерполяционный многочлен второй степени с узлами

Интегрируя этот многочлен, получаем элементарную формулу Симпсона:

Применяя эту формулу на каждом отрезке, получаем:

Остаточный член формулы Симпсона:

 

 

Решение.

Для варианта 4 функция .

 

С помощью остаточным членов интегрирования оценим шаги интегрирования, при которых величина погрешности каждой квадратурной формулы будет меньше

1. Формула трапеций:

 

 

 

 

 

2. Формула Симпсона:

 

 

 

 

 

С помощью программы вычисляем значения интеграла по квадратурным формулам с найденным шагом:

 

Результат формулы трапеций:

26.6713173333738

Результат формулы Симпсона:

26.6713173333867

 

Вычислим аналитически точное значение интеграла:

 

 

Найдем абсолютные погрешности результатов, полученных с помощью квадратурных формул:

 

Для формулы трапеций:

0,000000000040466= 4*10-11

Для формулы Симпсона:

0,000000000053366= 5.3*10-11

 

Измерим время вычисления интеграла по каждой из двух квадратурных формул ( миллисекундах). Для более наглядного значения произведем операцию вычисления 100000 раз. Полученные результаты:

 

Вывод: Формула Симпсона намного эффективнее формулы трапеций, т.к.  шаг метода Симпсона на 3 порядка больше, следовательно и время выполнения алгоритма метода Симпсона меньше.

 

Код программы на C#:

 

public task_8()

        {

            InitializeComponent();

            double h = 2 * Math.Pow(10, -6);

            lb1.Content += " Результат=" + Trapezoidal_rule(h , new Stopwatch()).ToString();

 

            h = 0.005;

 

            lb2.Content += " Результат=" + Simpson_rule(h, new Stopwatch()).ToString();

        }

 

 

 

 

        //метод трапеции

        public static double Trapezoidal_rule(double h, Stopwatch stopWatch)

        {

            double[] xmas = new double[(int)(0.8 / h)+1];

            double sum = 0;

            stopWatch.Start();

            for (int i = 0; i < xmas.Length; i++)

            {

                xmas[i] = 1 + h * i;

                if (i > 0 && i < xmas.Length - 1) sum += func(xmas[i]);

            }

 

 

            double trapec = h * ((func(xmas[0]) + func(xmas[xmas.Length - 1])) / 2 + sum);

            stopWatch.Stop();

            return trapec;

        }

 

        //подынтегральная функция

        public static double func(double x)

        {

            return 5.4 + 2.1 * x + 0.3 * Math.Pow(x, 2) + 2.1 * Math.Pow(x, 3) + 1.6 * (Math.Pow(x, 4) + Math.Pow(x, 5));

        }

 

        //метод симпсона

        public static double Simpson_rule(double h, Stopwatch stopWatch)

        {

            double[] xmas = new double[(int)(0.8 / h)+1];

            double sum = 0, sum1 = 0;

            stopWatch.Start();

            for (int i = 0; i < xmas.Length; i++)

            {

                xmas[i] = 1 + h * i;

                if (i < xmas.Length - 1)

                {

                    sum += func(xmas[i] + h / 2);

                    if (i > 0) sum1 += func(xmas[i]);

                }

            }

            double Simpson= h / 6 * (func(xmas[0]) + func(xmas[xmas.Length - 1]) + 4 * sum + 2 * sum1);

            stopWatch.Stop();

            return Simpson;

        }

 

        private void Button_Click(object sender, RoutedEventArgs e)

        {

 

            double h = 2 * Math.Pow(10, -6);

            Stopwatch stopWatch = new Stopwatch();

            for (int k = 0; k < 100; k++)

            {

 

                for (int i = 0; i < 1000; i++)

                {

                    Trapezoidal_rule(h, stopWatch);

                }

                ;

            }

            TimeSpan time = stopWatch.Elapsed;

            //lb1.Content += " Время вычисления=" + ((double)(count / 100000)).ToString();

            lb1.Content += " Время вычисления=" + ((double)(time.TotalMilliseconds / 100000)).ToString();

 

            h = 0.005;

            stopWatch = new Stopwatch();

            for (int i = 0; i < 100000; i++)

            {

                Simpson_rule(h,stopWatch);

            }

            time = stopWatch.Elapsed;

            lb2.Content += " Время вычисления=" + ((double)(time.TotalMilliseconds / 100000)).ToString();

        }

 

 

Задача 8.2. Исследовать поведение погрешности  приближенно  вычисленного  интеграла при уменьшении шага интегрирования.

 

Теоретический материал: Для погрешности квадратурной формулы при достаточно малом h справедливо приближенное равенство:

 

Уменьшение шага в 2 раза приводит к уменьшению погрешности примерно в 2k раз:

 

Вычитая из второго уравнения первое, приходим к приближенной формуле, называемой правилом Рунге:

 

 

Вычисляя оба значения: можно получить уточненное значение интеграла:

 

Для метода центральных прямоугольников k =2.

 

 

Решение.

Вычислим интеграл аналитически:

Требуется вычислить интеграл методом центральных прямоугольников.

Формула центральных прямоугольников: S= ; остаточный член    

Вычислим значения и построим по ним уточненное значение интеграла Iy.

Пример вычислений и графики погрешностей Iy в зависимости от шага h:

Код программы на C#:

        private double a = 0, b = 4;

        public Window1()

        {

            InitializeComponent();

            Bitmap bmp1 = new Bitmap(1000, 1000);

            Bitmap bmp2 = new Bitmap(1000, 1000);

            Graphics g1 = Graphics.FromImage(bmp1);

            Graphics g2 = Graphics.FromImage(bmp2);

            double h = (b - a) / 2.0;

            double I1, I2, Iy;

            double R1, R2, R3;

           

            double shir=(double)(2*(bmp1.Width-100)/(b-a));

            I1 = Central_Rectangle(h);

            h /= 2;

            I2 = Central_Rectangle(h);

            Iy = I2 + (I2 - I1) / 3;

            double vis =(double) ((bmp1.Height-100) / pogr(I2));

            System.Drawing.Point p1 = new System.Drawing.Point((int)(50+2*h * shir),bmp1.Height-50 - (int)(pogr(I2)*vis));

            System.Drawing.Point p3 = new System.Drawing.Point((int)(50 + 2*h * shir), bmp1.Height - 50 - (int)(pogr(Iy)*vis));

 

            System.Drawing.Pen pen = new System.Drawing.Pen(System.Drawing.Brushes.Red, 10);

            System.Drawing.Pen black_pen = new System.Drawing.Pen(System.Drawing.Brushes.Black, 3);

            System.Drawing.Pen black_pen_for = new System.Drawing.Pen(System.Drawing.Brushes.Black, 10);

            g1.DrawLine(black_pen_for,0,bmp1.Height - 50,bmp1.Width,bmp1.Height - 50);

            g2.DrawLine(black_pen_for,0,bmp2.Height - 50,bmp2.Width,bmp2.Height - 50);

            g1.DrawLine(black_pen_for,50,0,50,bmp1.Height);

            g2.DrawLine(black_pen_for,50,0,50,bmp2.Height);

 

            for (double x = 0; x <bmp1.Height-50; x+=10)

            {

                g1.DrawLine(black_pen, 0, bmp1.Height - 50 - (int)(x * vis), bmp1.Width, bmp1.Height - 50 - (int)(x * vis));

                g2.DrawLine(black_pen, 0, bmp2.Height - 50 - (int)(x * vis), bmp2.Width, bmp2.Height - 50 - (int)(x * vis));

            }

            for (double i = 0.1; i < 4; i+=0.1)

            {

                g1.DrawLine(black_pen, 50 + (int)(i * shir), 0, 50 + (int)(i * shir), bmp1.Height);

                g2.DrawLine(black_pen, 50 + (int)(i * shir), 0, 50 + (int)(i * shir), bmp2.Height);

            }

            for (; h > 0.0001; )

            {

                I1 = Central_Rectangle(h);

                h /= 2;

                I2 = Central_Rectangle(h);

                Iy = I2 + (I2 - I1) / 3;

                System.Drawing.Point p2=new System.Drawing.Point((int)(50+2*h*shir),bmp1.Height-50 -(int)(pogr(I2)*vis));

                System.Drawing.Point p4 = new System.Drawing.Point((int)(50+2*h * shir), bmp2.Height-50-(int)(pogr(Iy)*vis));

                g1.DrawLine(pen, p1,p2);

                g2.DrawLine(pen, p3,p4);

                p1 = p2;

                p3 = p4;

            }

            ImageSource imgSourceFromBitmap1 = System.Windows.Interop.Imaging.CreateBitmapSourceFromHBitmap((bmp1).GetHbitmap(),

                IntPtr.Zero, Int32Rect.Empty, BitmapSizeOptions.FromEmptyOptions());

            im1.Source = imgSourceFromBitmap1;

            ImageSource imgSourceFromBitmap2 = System.Windows.Interop.Imaging.CreateBitmapSourceFromHBitmap((bmp2).GetHbitmap(),

                IntPtr.Zero, Int32Rect.Empty, BitmapSizeOptions.FromEmptyOptions());

            im2.Source = imgSourceFromBitmap2;

        }

 

        public double Central_Rectangle(double h)

        {

            double[] xmas = new double[(int)((b-a) / h) + 1];

            double sum = 0;

            for (int i = 0; i < xmas.Length; i++)

            {

                xmas[i] = a + h * i;

                if ( i < xmas.Length - 1) sum += funct(xmas[i]+h/2);

            }

 

 

            return h * sum;

        }

 

        private double pogr (double Integral)

        {

            return Math.Abs( Integral + 512.5017427420596);

        }

        public double funct(double x)

        {

            return Math.Exp(2*x)*Math.Sin(x);

        }

 

        private void Button_Click(object sender, RoutedEventArgs e)

        {

            double h = Convert.ToDouble(tb.Text);

            double I1 = Central_Rectangle(h);

            double I2 = Central_Rectangle(h / 2);

            double Iy = I2 + (I2 - I1) / 3;

            r1.Content = "I(h)=" + I1 + " R1=" + pogr(I1);

            r2.Content = "I(h/2)=" + I2 + " R2=" + pogr(I2);

            r3.Content = "Iy=" + Iy + " R3=" + pogr(Iy);

        }

 

Задача 8.3. Вычислить интеграл от функции двух переменных с  точностью 0.01.

 

Теоретический материал.

Для интеграла от функции одной переменной формула левых прямоугольников выглядит так:

 

Формула правых прямоугольников:

 

 

 

Выведем формулу правых прямоугольников по x для элементарного прямоугольника .

Дано: интеграл   , N=11- номер варианта.

Методы вычисления: метод правых прямоугольников по x.

  метод левых прямоугольников по y.

 

 

 

 

 

Тест:

 

 

 

Исходя из элементарной формулы, строим составную, разделив исходный прямоугольник на n элементарных.

 

 

 

 

Решение:

Код программы на C#:

 

delegate void UpdateProgressBarDelegate(DependencyProperty dp, object value);

    public partial class task_8_3 : Window

    {

       

          //пределы интегрирования

        private double ax = 1, bx = 3, ay = 0, by = 2,hx,hy;

 

        // бэкграундный процесс

        BackgroundWorker worker;

 

        private decimal Integral1, Integral2;

        //флаг остановки вычисления

        bool stopprocess;

        public task_8_3()

        {

            InitializeComponent();

 

            worker = new BackgroundWorker();

            worker.WorkerSupportsCancellation = true;   //поддержка остановки процесса

            worker.DoWork += new DoWorkEventHandler(worker_DoWork); //добавление выполняемой процедуры

        }

 

        void worker_DoWork(object sender, DoWorkEventArgs e)

        {

             Integral1 = Method(hx, hy, pb1);

             Integral2 = Method(2 * hx, 2 * hy, pb2);  

 

        }

 

        //подынтегральная функция

        private decimal func(double x,double y)

        {

            return (decimal)(Math.Pow(11, x + y) * Math.Sin(x * y));

        }

 

        private void prbarUpdate(ProgressBar a)

        {

            a.Value += 1;

        }

 

        // по x - правых прямоугольников, по y - левых

        public decimal Method(double hx, double hy, ProgressBar a)

        {

            double[] xmas = new double[(int)((bx - ax) / hx) + 1];

            double[] ymas = new double[(int)((by - ay) / hy) + 1];

            //создаем делегат для обновления прогресбара

            UpdateProgressBarDelegate updProgress = new UpdateProgressBarDelegate(a.SetValue);

 

            //инициализация массивов x,y

            for (int i = 0; i < xmas.Length; i++)

            {

                xmas[i] = ax + hx * i;

            }

            for (int i = 0; i <ymas.Length; i++)

            {

                ymas[i] = ay + hy * i;

            }

 

            decimal sum = 0;

 

            double value = 0;

            for (int i = 1; i < xmas.Length && !stopprocess; ++i)

            {

                //сумма всех значений подынтегральной  функции при фиксированном x(i) для всех y(j), кроме последнего

                sum += ParallelEnumerable.Range(0, ymas.Length - 1).Sum(j => func(xmas[i], ymas[j]));

 

                //обновляем прогресбар

                Dispatcher.Invoke(updProgress, new object[] { ProgressBar.ValueProperty, ++value });

            }

 

 

            return sum * ((decimal)(hx * hy));

        }

 

        private void Button_Click(object sender, RoutedEventArgs e)

        {

            lb1.Content = null;

            lb2.Content = null;

 

            stopprocess = false;

 

            hx = Convert.ToDouble(tbx.Text);

            hy = Convert.ToDouble(tby.Text);

 

            //изменяем максимальные значения прогрессбаров

            pb1.Maximum = (int)((bx - ax) / hx) + 1;

            pb2.Maximum = (int)((bx - ax) / (2 * hx)) + 1;

 

            //обнуляем прогрессбары

            pb1.Value = 0;