Информатика_1 / C / lecture7 / lecture7

Лекция 7. Численные методы вычисления определенных интегралов Table of contents. 1. Метод прямоугольников для вычисления интегралов. 2. Многочлены Лагранжа и метод Симпсона для вычисления интегралов. 3. Метод Монте-Карло для вычисления интегралов. 4. Задачи. Численные методы вычисления определенных интегралов. В этой лекции рассматриваются численные методы для вычисления определенных интегралов. Будут рассмотрены методы прямоугольников, Симпсона и Монте-Карло. Метод прямоугольников для вычисления определеннных интегралов . В методе прямоугольников для численного вычисления интеграла используется геометрическая интерпретация определенного интеграла , как площади под кривой. Для вычисления этой площади весь интервал интегрирования [a,b] разбивается на n равных подинтервалов длины h=(b-a)/n. Площадь под подинтегральной кривой приближенно заменяется на сумму площадей прямоугольников, как это показано на следующем рисунке. Сумма площадей всех прямоугльников вычисляется по формуле Разница между точным значением интеграла I и приближенным , представленным интегральной суммой S зависит от шага интегрирования h cледующим образом Из этой формулы видно, что погрешность вычисления интеграла определяется величиной шага интегрирования h. Чем меньше шаг интегрирования , тем точнее интегральная сумма S аппроксимирует значение интеграла I. Исходя из этого строится алгоритм для вычисления итеграла с заданной точностью . Считается , что интегральная сумма S представляет значение интеграла I c точностью eps если разница по абсолютной величине между интегральными суммами S0 и S подсчитенными для h и h/2 cсоответственно не превышает eps. Ниже приведена блок-схема для вычисления определенного интеграла с заданной точностью eps. В этой блок-схеме при числе разбиений интервала интегрирования n большем чем nmax считается , что сходимость не достигнута и делаетя аварийный выход. Далее приводится пример программы , в которой функция rect вычисляет интеграл по методу прямоугольников и пример ее вызова из функции main (файл l7_1.c): /*1*/ #include<stdio.h> /*2*/ #include<math.h> /*3*/ int rect(double,double,double,double (*)(double),double *); /*4*/ double func(double); /*5*/ void main(void){ /*6*/ double s,eps=0.001,a=0,b=1; /*7*/ int error; /*8*/ error=rect(a,b,eps,func,&s); /*9*/ if(error==0)printf("integral=%e\n",s); /*10*/ else printf("error in rect\n"); /*11*/ } /*12*/ int rect(double x1,double x2,double e,double (*g)(double),double *s) /*13*/ {int n=20,nmax=16384,i; double h,s0; /*14*/ h=(x2-x1)/n; /*15*/ for(i=1,*s=0;i<=n;i++)*s+=g(x1+(i-1)*h); /*16*/ *s*=h; /*17*/ do { s0=*s; n*=2; h/=2; /*18*/ for(i=1,*s=0;i<=n;i++)*s+=g(x1+(i-1)*h); /*19*/ *s*=h; /*20*/ if(fabs(*s-s0)<e)return 0; /*21*/ } while(n<nmax); /*22*/ //аварийный выход /*23*/ return 1; /*24*/ } /*25*/ double func(double t){ return t*t; } В строках 12-24 описана функция rect , предназначенная для вычисления определенных интегралов методом прямоугольников по блок-схеме , приведенной выше. Формальные аргументы: x1,x2 - нижний и верхний пределы интегрирования, e - погрешность вычисления интеграла, g - имя подинтегральной функции, s - адрес вещественной переменнной для результата вычисления интеграла. С помощью оператора return функция rect возвращает код ошибки, 0 - интеграл вычислен, 1 - число разбиений интервала интегрирования превысило максимально допустимое значение (nmax), интеграл не вычислен с требуемой точностью. В строке 8 вызывается функция rect , код возврата записывается в переменную error. В случае успешного завершения функции rect в переменной error хранится ноль и тогда значение интеграла берется из переменной s и выводится на консоль ( строка 9), в противном случае на консоль выдается сообщение об ошибке. Интерполяционные многочлены Лагранжа и метод Симпсона для вычисления определенных интегралов . Для вывода формулы Симпсона для интегральной суммы введем определение интерполяционного многочлена Лагранжа для функции f(x). Интерполяционным многочленом степени n для функции f(x) называется такой многочлен степени n который удовлетворяет следующим условиям : Один из вариантов интерполяционного многочлена степени n для функции f(x) был дан Лагранжем В методе Симпсона для вычисления определенного интеграла весь интервал интегрирования [a,b] разбивается на четное число подинтервалов равной длины h=(b-a)/n. Затем на каждой паре соседних подинтервалов подинтегральная функция f(x) заменяется многочленом Лагранжа второй степени, т. е. Ниже приводится пример функции , вычисляющей определенный интеграл по методу Симпсона. Функция simpson работает по тому-же алгоритму (блок-схеме), что и функция rect из предыдущего параграфа, но интегральная сумма вычисляется по формуле Симпсона 7.5. (код функции хранится в файле l7_2.c) /*1*/ int simpson(double a,double b,double e,double (*g)(double),double *s) /*2*/ { int n=20,nmax=16384,i; double h, s0 ; /*3*/ h=(b-a)/n; /*4*/ for( i=1, *s=0 ; i<=n-1 ; i++ ) if ( i % 2 == 0 ) *s += 2*g( a+i*h) ; else *s+=4*g( a+i*h ); /*5*/ *s+=g(a)+g(b); *s*=h/3; /*6*/ do { s0=*s; n*=2; h/=2; /*7*/ for(i=1,*s=0;i<=n-1;i++)if(i%2==0)*s+=2*g(a+i*h); else *s+=4*g(a+i*h); /*8*/ *s+=g(a)+g(b); *s*=h/3; /*9*/ if(fabs(*s-s0)<e)return 0; /*10*/ } while(n<nmax); /*11*/ //аварийный выход /*12*/ return 1; /*13*/ } Формальные параметры функции: a,b - пределы интегрирования, e - допустимая погрешность при вычислении интеграла, g - подинтегральная функция, s - адрес (указатель) вещественной переменнной для ычисленного значения интеграла. Функция возвращает с помощью оператора return код завершения , 0 - если интеграл вычислен с заданной точностью и 1 - если число разбиений интервала интегрирования превысило максимально допустимое значение nmax. Функция simpson отличается от функции rect из предыдущего параграфа лишь заложенной формулой для интегральной суммы. В данном случае формула Симпсона для интегральной суммы заложена в строках 4,5 и 7,8. Отмечу два момента , связанных с приведенным примером. Во-первых , последним аргументом функции simpson является адрес вещественной переменной типа double, предназначенный для записи результата вычисления интеграла. Поскольку s указатель, то для получения значения , хранящегося по этому адресу применяется операция разадресации *. И во-вторых, при вычислении интегральной суммы 7.5 в строках 4 и 7 применяется операция получения остатка от целого деления %. C помощью этой опреации при четных i в интегральную сумму *s суммируется значение подинтегральной функции с множителем 2, а при нечетных i - с множителем -4. Слушателям в качестве упражнения рекомендуется написать вызывающую программу для функции simpson. Метод Монте-Карло для вычисления определенных интегралов . Метод Монте-Карло занимает особое положение среди методов вычисления определенных интегралов по двум причинам. Во-первых, это единственный метод позволяющий вычислять интегралы высокой кратности. И во-вторых, это метод, который дает лишь вероятностные гарантии степени точности вычисления интегралов. В этой лекции приводится описание простейшего варианта метода и без математических доказательств. Итак, для вычисления однократного интеграла методом Монте-Карло может быть применена формула Все что нужно для вычисления интегральной суммы по формуле 7.6 - зто научиться получать случайные числа равномерно распределенные на интервале [a,b]. Для этой цели можно использовать генератор случайных чисел , входящий в состав стандартных библиотек , поставляемых с компилятором. В заголовочном файле stdlib.h приведен заголовок функции int rand(void), генерирующей случайное целое число из интервала [0, RAND_MAX], где RAND_MAX - макрос, определенный в том-же заголовочном файле. С помощью функции rand и макроса RAND_MAX легко получить случайное вещественное число , равномерно распределенное на интервале [0,1] - например, результатом выполнения оператора (double)rand()/RAND_MAX является случайное вещественное число из интервала [0,1]. Опрерация (double) - приведения к типу double результата , возвращенного функцией rand в приведенном выше выражении необходима - иначе (в случае отсутствия этого приведения) будет выполняться целое деление т.к. и макрос RAND_MAX и функция rand принимают целочисленные значения , а в языке Си целое деление осуществляется нацело. Имея случайное вещественное число из интервала [0,1] легко получить случайное число из любого интервала. Например, если z - cлучайное число из инетервала [0,1] , тогда x=a+(b-a)*z - cлучайное число из интервала [a,b]. Как видно из приведенных выше оценок погрешности формулы 7.6 точность вычисления интеграла и в методе Монте-Карло определяется числом слагаемых N в интегральной сумме - чем больше слагаемых тем точнее результат . Замечу, что формула для интегральной суммы 7.6 отличается от аналогичной формулы метода прямоугольников лишь тем, что в методе прямоугольников мы суммируем значения подинтегральной функции в равноотстоящих точках x , а в методе Monte-Karlo - в случайных. Для вычисления определенного интеграла методом Монте-Карло с заданной точностью может быть применена та-же самая блок-схема, что применялась нами для вычисления интегралов методами прямоугольников и Симпсона. Ниже приведен пример функции с именем monte , вычисляющей определенный интеграл методом Монте-Карло . Список формальных параметров этой функции и возвращаемое значение совершенно аналогичны соответсвенным параметрам функций rect и simpson и поэтому не обсуждаются. int monte(double x1,double x2,double e,double (*g)(double),double *s) {int n=20,nmax=16384,i; double h,s0,x; h=(x2-x1)/n; for(i=1,*s=0;i<=n;i++){x=x1+(x2-x1)*(double)rand()/RAND_MAX;*s+=g(x);} *s*=h; do { s0=*s; n*=2; h/=2; for(i=1,*s=0;i<=n;i++){x=x1+(x2-x1)*(double)rand()/RAND_MAX;*s+=g(x);} *s*=h; if(fabs(*s-s0)<e)return 0; } while(n<nmax); //аварийный выход return 1; } Отмечу лишь, что для использования этой функции необходимо включить в программу заголовочный файл stdlib.h, в котором объявлен заголовок функции rand и макрос RAND_MAX, определяющий максимальное случайное число , генерируемое функцией rand. Метод Монте-Карло легко обобщается на интегралы произвольной кратности. Например двукратный интеграл может быть вычислен по формуле Ниже приводится пример функции с именем monte2 для вычисления двукратного интеграла int monte2(double a,double b,double c,double d,double e,double (*g)(double,double),double *s) {int n=20,nmax=16384,i; double h,s0,x,y; h=(b-a)*(d-c)/n; for(i=1,*s=0;i<=n;i++) {x=a+(b-a)*(double)rand()/RAND_MAX; y=c+(d-c)*(double)rand()/RAND_MAX; *s+=g(x,y);} *s*=h; do { s0=*s; n*=2; h/=2; for(i=1,*s=0;i<=n;i++){ x=a+(b-a)*(double)rand()/RAND_MAX; y=c+(d-c)*(double)rand()/RAND_MAX; *s+=g(x,y);} *s*=h; if(fabs(*s-s0)<e)return 0; } while(n<nmax); //аварийный выход return 1; } В списке формальных аргументов этой функции , по сранению с функцией monte , добавились два входных параметра - c и d верхний и нижний пределы второго интеграла. Кроме того подинтегральная функция g теперь описана , как функция двух переменных. Задачи к лекции . Задача 1. Написать функцию для вычисления интеграла методом Симпсона. Функция должна возвращать код ошибки и результат вычисления интеграла через список аргументов. Написать вызывающую функцию. Задача 2. Написать функцию для вычисления однократного интеграла методом Монте-Карло. Функция должна возвращать код ошибки и результат вычисления интеграла через список аргументов. Написать вызывающую функцию. Задача 3. Написать функцию для вычисления интеграла методом Симпсона. Функция должна возвращать код ошибки и результат с помощью оператора return, используя структуру. Написать вызывающую функцию. Задача 4. Обобщив метод Монте-Карло для случя трехкратного интеграла написать функцию для его вычисления. Код ошибки и значение интеграла должны возвращаться через список аргументов. Написать вызывающую прграмму. previous next home