Закрыть

Интегрирование в Maple.

Автор: Васин Алексей
Опубликовано 08.02.2008 в 15:56
Раздел: Прикладные пакеты

Краткий обзор

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

Если формальное дифференцирование было реализовано в 1953 году, то первые шаги в области формального интегрирования были сделаны только в 1961 году [см. Davenport 1987]. Чем обусловлена эта задержка?

Причину следует искать в существенном различии между этими операциями. В отличии от алгоритмической процедуры дифференцирования, интегрирование выглядит случайным набором правил. Кажется, что существует только одно общее правило:

"Интеграл от суммы функций равен сумме интегралов слагаемых" (да и это правило не всегда справедливо, если мы фиксируем пространства, которым должны принадлежать функции). Для операций, отличных от сложения и вычитания, общих правил не существует. Например, из того, что мы знаем, как интегрировать

x2 и ex, вовсе не следует, что мы знаем, как интегрировать их композицию. Таким образом, обычно все мы изучаем несколько "методов", таких, как интегрирование по частям, интегрирование с помощью подстановки и т.д. Кроме того, априори мы не знаем, какой метод или комбинация методов сработает для данного интеграла. Так что первые шаги основывались на той же эвристике, которую использует человек. Однако, этот способ был весьма скоро вытеснен (1967-1971) истинно алгоритмическими методами. К настоящему времени существует развитая теория интегрирования, и мы можем дать только совсем краткий ее обзор.

Практически наилучший алгоритм интегрирования на 1997 год реализован в системе Macsyma, решающей 95% примеров из справочника Камке.

Поскольку дифференцирование, несомненно, проще, чем интегрирование, имеет смысл переформулировать проблему интегрирования как "обратную задачу" для дифференцирования, т.е. для данной функции f вместо отыскания ее интеграла g мы ищем функцию g, такую, что g'=f. Определение: Для двух данных классов функций F и G задача интегрирования из F в G состоит в том, чтобы найти алгоритм, который для любого элемента f из F либо выдает элемент g из G, такой, что f=g', либо доказывает, что в G не существует такого элемента g, что f=g'.

 

В 1968 году Ричардсон доказал теорему, утверждающую, что проблема интегрирования для F=G=Q неразрешима, где Q поле рациональных функций от функций i, Pi , exp , log , и ||-абсолютная величина. Но это не должно разочаровывать, в действительности в этом поле не может быть решена уже проблема определения, является ли данная константа нулем (и, как следствие, не может быть реализован даже алгоритм деления). Реально, вводят понятия неопределимой константы, поля эффективных констант и множество других достаточно сложных конструкций.

Только в 1989 году удалось доказать, что интегрирование является алгоритмической процедурой для весьма ограниченного класса функций. Однако, остается еще задача о построении достаточно эффективных алгоритмов интегрирования. Фактически эта задача остается на переднем крае и математики и возможностей систем компьютерной алгебры.

Задачу интегрирования можно рассматривать как решение простейшего дифференциального уравнения diff(y,x)=f. Замечания, сделанные выше об алгоритме интегрирования остаются справедливыми и для алгоритма решения дифференциальных уравнений. Например, рассмотрим уравнение y' + f y = g. Хорошо известен метод решения уравнения данного типа - надо воспользоваться подстановкой y = z exp ( - Int ( f,x) ).

Получаем g = y' + f y = z' exp( - Int f ) таким образом решение равно y = exp ( -Int f ) Int g exp ( Int f ). Здесь Int обозначает интеграл по х.

В общем случае этот метод решения алгоритмически неудовлетворителен, поскольку алгоритм интегрирования переформулирует этот интеграл опять в виде уравнения, с которого мы начинали.

В силу этого, задача опять формулируется следующим образом: для двух данных функций f и g (какого-то фиксированного класса) найти функцию y, такую что y' + f y = g или доказать, что такой функции в данном классе нет.

Современные алгоритмы решения дифференциальных уравнений достаточно сложны, и, вероятно, наиболее полно реализованны в системе Macsyma, Maple и Mathematica.

Решения дифференциального уравнения в конечном виде, если они существуют, очень полезны. Но у большинства дифференциальных уравнений в физике таких решений нет. Именно поэтому в системах символьных вычислений реализованны алгоритмы не только численного решения уравнений, но и алгоритмы решений в виде рядов (например, асимптотические решения). Математическое обоснование существующих алгоритмов является целой отдельной областью математики, которую мы не будем затрагивать в данной статье.

 

Встроенные команды интегрирования Maple.

В стандартной (загруженной по умолчанию) библиотеке находяться процедуры

Int ( expr, options ) и int ( expr,options ), т.е. символ и значение.

Синтаксис и опции можно найти в справочнике int.

 

Пример 1 (определенные и неопределенные интегралы)

f:=sin(x):

g:=Int(f,x);

value(g);

int(f,x);

 

Таким образом, мы видим. что обозначения максимально приближены к обычной математической записи. Для вычисления определенного интеграла надо дополнительно указать пределы интегрирования:

 

f:=x*exp(x):

g:=Int(f,x=0..1);

value(g);

int(f,x=0..1);

Рассмотрим более сложные интегралы

f:=sin(x^3):

g:=int(f,x=0..Pi);

h:=evalf(g);

или

g:=Int(x/(x^3-x^2+1),x);

value(g);

 

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

 

restart: readlib(branches):

branches(ln);

 

Рассмотрим интегралы с параметрами

 

f:=exp(-a*x)*ln(x)*sqrt(x);

int(f,x=0..infinity);

 

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

 

assume(a<0): int(f,x=0..infinity);

Сравним с

assume(a>0): int(f,x=0..infinity);

 

Полученные при интегрировании эллиптические функции, зависящие от радикалов можно упростить

 

answer := int( 1/sqrt( sin(x) ), x=0..Pi/2 );

radnormal(answer);

 

Интегралы, которые не выражаются через алгебраические функции, можно представить как ряд

 

int(exp(x^3), x );

series(%, x=0);

 

В общем случае синтаксис для численного интегрирования имеет вид

 

evalf ( Int ( f, x = a..b , digits , flag ) ),

 

где digits- точность вычисления, и flag- код численного метода, см. int[numerical]

 

int( exp(v-v^2/2)/(1+1/2*exp(v)), v = -infinity..infinity );

evalf(", 20);

 

Здесь digits равно 20. Другой пример - интегрирование гамма функции с использованием алгоритма Ньютона-Котеса:

 

evalf(Int(1/GAMMA(x), x = 0..2, 10, _NCrule));

 

Более сложные опции позволяют вычислять главное значение интеграла и реализуют некоторые дополнительные возможности

 

int(1/x^3, x=-1..2, 'CauchyPrincipalValue');

int( 1/(x+a)^2, x=0..2, 'continuous' );

 

В принципе, различные эвристические методы, такие как интегрирование по частям и замены переменных, не используются явно, например

 

Пример 2 (интегрирование по частям и другие методы)

 

Рассмотрим несколько интегралов, которые берутся по частям или с использованием различных замен переменных

 

int(x*ln(x),x);

int(x*sin(x^2),x);

Int(x^2/(sqrt(1-x^2))^(3/2),x):

value(");

 

Тем не менее, для обучения, мы можем использовать подобные команды явно. Эти команды содержатся в пакете student

 

Пример 3 (интегрирование с пакетом student)

 

Для начала мы очистим рабочую память и загрузим данный пакет.

 

restart:

with(student);

 

При загрузке любого пакета Maple выдает список всех загружаемых команд.

 

Замены переменных:

f:=Int(sqrt(1-x^2),x=a..b);

g:=changevar(x=sin(u),f,u);

 

Здесь мы вводим новую независимую переменную интегрирования u, связанную со старой переменной уравнением x=sin(u) (пределы интегрирования изменились автоматически!) Посмотрим на результат

 

value(g);

 

В двойных интегралах интегралах это позволяет, например, перейти к полярным координатам:

 

changevar({x=r*cos(t),y=r*sin(t)}, Doubleint(1,x,y),[t,r] );

 

Интегрирование по частям:

напомним, что правило интегрирования по частям записывается следующим образом

 

f:=Int (u,v)=u*v-Int(v,u);

 

Синтаксис данной команды intparts(f,u)

 

a:=Int(x^2*exp(x),x);

 

Проинтегрируем по частям

 

b:=intparts(a,x^2);

 

проинтегрируем по частям второй раз

 

c:=intparts(b,2*x);

 

и получим результат

 

d:=value(c);

 

Другие полезные команды для интегрирования из данного пакета

 

integrand(Int(h(x),x)); и

Doubleint(h(x,y),x,y);

 

Tripleint(h,x=1..n,y=2..4,z=w..u);

 

Кроме этого, в пакете student реализовано несколько численных методов: leftbox, leftsum , rightbox, rightsum, middlebox, middlesum simpson, trapezoid.

Синтаксис для всех команд одинаков, рассмотрим только один пример

 

Пример 4 (численное интегрирование)

 

restart: with(student):

f:=x^5*ln(x);

 

Явное интегрирование данной функции можно посмотреть в справочнике, мы же сравним численные значения

 

ex:=evalf(int(f,x=1..3));

s:=evalf(simpson(f,x=1..3,16));

t:=evalf(trapezoid(f,x=1..3,16));

ls:=evalf(leftsum(f,x=1..3,16));

 

Покажем это на рисунке

 

leftbox(f,x=1..3,16,color=red,shading=BLUE);

 

Аналогично и для других команд

 

rs:=evalf(rightsum(f,x=1..3,16));

 

rightbox(f,x=1..3,16,color=red,shading=BLUE); и

ms:=evalf(middletsum(f,x=1..3,16)));

middlebox(f,x=1..3,16,color=red,shading=BLUE);

 

Как говорил Козьма Прутков:-"Нельзя объять необятное". В силу этого, тяжело загруженный универсальностью грузовик Maple иногда довольно долго и мучительно пробирается по стандартным наезженым дорогам математики. В этом случае полезно пожертвовать универсальностью, вспомнить что Maple является языком программирования и написать свою собственную программу. Приведем пример:

 

Пример 5 (численное интегрирование - метод Гаусса)

 

Напомним, что интегрирование с помощью замены переменной приводится к интегрированию на отрезке [0,1], а затем используется квадратурная формула Гаусса. Мы зададим руками соответствующие абциссы и веса (8 или 16 точек) и перепишем стандартную программу в виде

 

restart:

dg2:=proc(f,a,b,tol)

local X,W,co,aa,bb,c1,c2,s8,s16,i,u,res,ndig;

W:=[.1012285362903763,

.2223810344533745,

.3137066458778873,

.3626837833783619,

.02715245941175409,

.06225352393864789,

.09515851168249278,

.1246289712555339,

.1495959888165767,

.1691565193950025,

.1826034150449236,

.1894506104550685]:

X:=[.9602898564975362,

.7966664774136267,

.5255324099163289,

.1834346424956498,

.9894009349916499,

.9445750230732326,

.8656312023878317,

.7554044083550030,

.6178762444026437,

.4580167776572274,

.2816035507792589,

.09501250983763744]:

if b=a then RETURN(0) fi;

res:=0: ndig:=16; #Заметим если ndig>16 требуются совсем другие X,W!

co:=evalf(0.005/(b-a),ndig);

bb:=a;

while evalf(bb,ndig)<evalf(b,ndig) do

aa:=bb;

bb:=b;

while true do

c1:=(bb+aa)/2;

c2:=(bb-aa)/2;

if evalf(1+abs(co*c2),ndig) = evalf(1,ndig) then

RETURN(`Exceeded accuracy limit`,res);

fi;

s8:=0;

for i from 1 to 4 do

u:=evalf(c2*X[i],ndig);

s8:=evalf(s8+W[i]*(f(c1+u)+f(c1-u)),ndig);

od:

s8:=evalf(c2*s8,ndig);

s16:=0;

for i from 5 to 12 do

u:=evalf(c2*X[i],ndig);

s16:=evalf(s16+W[i]*(f(c1+u)+f(c1-u)),ndig);

od:

s16:=evalf(c2*s16,ndig);

if evalf(abs(s16-s8),ndig) < evalf(tol*(1+abs(s16)),ndig)

then break fi;

bb:=c1;

od:

res:=evalf(res+s16,ndig);

od:

res;

end:

 

Мы не будем вдаваться в подробности численного интегрирования. На данном этапе нашей целью является сравнение скорости работы данной программы и встроеной универсальной Maple-программы. Для примера проинтегрируем функцию 

 

fun:=x->x^31;

 

В переведенной нами программе последний аргумент задает точность вычислений

 

dg2(fun,0,5,0.000000000001);

 

Теперь сравним с точным значением

 

evalf(%-(5^32)/32,16); и посмотрим как быстро произведет вычисления Maple-программа. Для оценки времени выполнения программ воспользуемся встроенной функцией Maple: time(e), где е - некоторое выражение или функция.

 

time(evalf(Int(fun(x),x=0..5),20));

time(dg2(sin,0,19*Pi,0.001));

 

Если у Вас не суперкомпьютер, то Вы сами убедитесь в медлитольности универсальной Maple-программы. На моем собственном компьютере приведенная выше программа  работает почти в 6 раз быстрее стандартной программы Maple.

 

Посмотрим еще несколько примеров и убедимся как от задаваемой точности зависит скорость работы и приближение к точному значению

 

f:=x->1/x;

dg2(f,1,50,0.1);

evalf(%-ln(50),16);

 

dg2(f,1,50,0.001);

evalf(%-ln(50),16);

 

dg2(f,1,50,0.0001);

evalf(%-ln(50),16);

 

dg2(f,1,50,0.000001);

evalf(%-ln(50),16);

 

dg2(x->1/x,1,50,0.0000000000000001);

evalf("-ln(50),16);

 

Двойные и тройные интегралы в Maple. Для примера рассмотрим несколько задач с первого курса.

 

Пример 6 (метод Монте-Карло для вычисления определенных интегралов).

 

В данном примере мы рассмотрим на простейшем примере основную идею метода Монте-Карло.

Рассмотрим единичную окружность внутри квадрата с вершинами (-1,-1), (1,-1), (1,1), (-1,1)

 

restart:

with(plots):

P1 := polygonplot([[-1,-1],[1,-1],[1,1],[-1,1]],color=green):

P2 := polarplot(1,color=gold):

display([P1,P2],scaling=constrained);

 

Предположим, что это мишень в которую мы стреляем. Некоторые наши стрелы попадут внутрь круга, некоторые остануться вне круга. Далее мы можем преположить, что отношение числа стрел попавших внутрь круга NR к общему числу выпущенных нами стрел N будет пропорционально отношению площадей круга и квадрата. Так как площадь квадрата нам известна S=4, то можно оценить площадь круга как S=4*NR/N.

Грубо говоря, в приведенном выше рассуждении содержиться основная идея метода Монте-Карло.

Посмотрим теперь как эту идею можно воплотить с помощью Maple. Итак, мы стреляем и попадаем в точку с координатами (x,y) которые создает генератор случайных чисел

 

x := evalf(rand(10001)/5000-1):

y := evalf(rand(10001)/5000-1):

 

После тысячи выстрелов

 

N := 1000;

NR := 0;

 

подсчитаем количество попаданий внутрь круга

 

for i from 1 to N do

if x()^2 + y()^2 < 1 then NR := NR +1 fi

od:

 

и вычислим приближение к площади единичного круга, равной числу

 

pi=3,14;

4*evalf(NR/N);

 

Конечно, увеличивая число выстрелов мы получим более точное приближение.

Найти объем тела ограниченного поверхностями

 

f=x^2+2y^2 - снизу, и

g=6-2x^2-y^2 - сверху,

 

единичным квадратом

 

0 <= x ,

x<= 1,

0 <= y ,

y <= 1.

 

Нарисуем поверхности

restart:

with(plots):

 

f:=x^2/8+y^2/4: g:=1-x^2/4-y^2/8:

surf:=plot3d({ f,g },x=0..1.2,y=0..1.2,lightmodel=light1):

planes:=plot3d( {[0,u,z],[1,u,z],[u,1,z]},u=0..1,z=0..1, grid=[20,20],orientation=[-65,67],style=patch,lightmodel=light1,axes=boxed):

 

Выведем их на экран

 

display3d({surf,planes});

 

Полный объем единичного куба равен V=1. Генератором случайных чисел зададим координаты внутри куба

 

X :=proc() evalf(rand(10000)/9999): end;

Y :=proc() evalf(rand(10000)/9999): end;

Z :=proc() evalf(rand(10000)/9999): end;

 

Теперь сделам тысячу выстрелов

 

N := 1000:

NR := 0:

 

и подсчитаем число попаданий

 

for i from 1 to N do

if Z() > (X()^2/8+Y()^2/4) then

if Z() < (1-X()^2/4-Y()^2/8) then NR := NR +1 fi:

fi:

font color="#ff0000">od:

 

Итого, объем тела примерно равен

 

evalf(NR/N);

 

Теперь сравним с точным значением

 

Int(Int(g-f,y=0..1),x=0..1):

evalf(value(%));

 

Примеры решения задач на интегрирование

 

Задача 1. Найти площадь ограниченную параболой x=2y^2 и линией x=y+3.

 

Решение 1.

Нарисуем эту область

 

restart:

x1:=2*y^2: x2:=y+3:

plot([x1,x2],y=-2..2);

 

Найдем точки пересечения кривых (пределы интегрирования)

 

solve({x=2*y^2,x=y+3});

 

s:=Int(Int(1,x=2*y^2..y+3),y=-1..3/2);

value(s);

 

Задача 2. Найти среднее значение функции f=y ch (x) над областью, ограниченной кривыми x=0 , y=2 , y=sqrt(x).

 

Решение 2.

Посмотрим на эту область

 

restart:

plot({[0,t,t=-1..2.5], [t,2,t= 0..5], [t,sqrt(t),t=0..5]},thickness=2);

 

Область описывается, как 0<y<2, 0<x<y^2. Вычислим среднее

 

s:=Int(Int(y*cosh(x),x=0..y^2),y=0..2)/Int(Int(1,x=0..y^2),y=0..2);

 

value(s); evalf(%);

 

Задача 3. Найти объем тела ограниченного поверхностями f=x^2+2y^2 - сверху, и g=6-2x^2-y^2 - снизу, над кругом единичного радиуса x^2+y^2 <1.

 

Решение 3.

Нам необходимо найти объем следующей фигуры.

 

restart: with(plots):

 

f:=x^2+2*y^2:

g:=6-2*x^2-y^2:

 

Нарисуем поверхности

 

surf:=plot3d({ f, g },x=-1.2..1.2,y=-1.2..1.2):

 

и цилиндр единичного радиуса

 

cyl:=plot3d([cos(t),sin(t),z],t=0..2*Pi,z=0..6):

 

Выведем их на экран

 

display3d({surf,cyl});

 

Чтобы найти объем данной фигуры необходимо вычислить интеграл

 

Int(Int(g-f,y=-sqrt(1-x^2)..sqrt(1-x^2)),x=-1..1);

value(%);

 

Задача 4. Используя сферические координаты найти объем тела, ограниченного сферой с x^2+y^2+z^2=4вырезанным в ней цилиндром x^2+y^2=1.

 

Решение 4.

 

В сферических координатах сфера задается как r=2, цилиндр как r sin(f)=1 . Нам надо выразить декартовы координаты x,y,z через сферические, подставив данные значения r для сферы и для цилиндра

 

restart: with(plots):

rs:=2; rc:=1/sin(phi);

x:=sin(phi)*cos(theta);

y:=sin(phi)*sin(theta);

z:=cos(phi);

 

Зададим точку на сфере и цилиндре

 

sf:=[rs*x,rs*y,rs*z]:

cyl:=[rc*x,rc*y,2*phi-Pi]:

 

Нарисуем сферу и цилиндр

 

plot3d({sf,cyl},theta=0..2*Pi,phi=0..Pi);

 

Вычислим пределы интегрирования, решив систему уравнений

 

_EnvAllSolutions := true:

solve({r=2,r*sin(phi)=1});

 

здесь B1 и Z1 целые числа. Таким образом, наша фигура может быть получена вращением следующего сегмента вокруг оси Z

 

plot({[rs*sin(v),rs*cos(v),v=Pi/6..5*Pi/6], [1,v,v=-sqrt(3)..sqrt(3)]},-3..3,-2..2);

 

Вычислим ее объем

 

Int(Int(Int((rho^2)*sin(phi),

rho=1/sin(phi)..2),phi=Pi/6..5*Pi/6), theta=0..2*Pi);

value(%);

Комментарии (0)

Комментировать могут только зарегистрированные пользователи

Подразделы
Новые статьи
Aрхив статей