Математика. Настоящий учебник посвящен системе Mathematica прикладному пакету компьютерной алгебры, при помощи которого можно решать любые задачи, в которых в той или иной форме встречается математика
 Скачать 4.43 Mb.
|
|
is the answer, what is the question? — а в числителе (2n − 1)!/2 n −1 (n − 1)!? Иными словами, надеемся ли мы 115 увидеть что-то в духе 2n Y i=1 i ( −1) i −1 = (2n − 1)! 2 2n −1 n!(n − 1)! ? Да нет же, мы уже говорили, что проводимые системами компьютерной алгебры вычисления основаны на совершенно других принципах, а реаль- ность превосходит все наши представления о ней. Если система компьютерной алгебры и сможет найти формулу для этого произведения, то это, конечно, будет совсем другая формула! Либо та же формула, но найденная как-то иначе и записанная в совершенно другом виде. Абсо- лютно принципиальный момент, который сразу очень многое объясняет в истории математики, состоит в том, что Mathematica скорее найдет со- вершенно невероятную формулу в духе Эйлера и Рамануджана, чем ту простую формулу, которую написал бы современный ма- тематик. Ясно, что никакой мистики здесь нет, она действует тем спосо- бом, который с ее точки зрения является наиболее экономным, разумным и простым! Это всего лишь значит, что ее представление об экономии мысли отличается от нашего!!! В любом случае, ценой небанального напряжения мысли (больше секун- ды раздумий!) Mathematica находит общую формулу. Для полного реализ- ма возьмем произведение до n, а не до 2n, чтобы не оставалось уже никаких сомнений, что начиная вычисление система не может знать ответа: In[121]:= Product[i^((-1)^(i-1)), {i,1,n}] Out[121]=E^((-1/2)*(-Log[2]+(-1)^n*Log[2]+Log[Pi])+ (-1)^n*(Log[Gamma[(1+n)/2]]-Log[Gamma[(2+n)/2]])) Ну как, все еще из Таблицы? Да, но почему тогда у Maple нет доступа к этой Таблице? Взяв теперь произведение до 2n, мы придем к формуле, которая хотя и эквивалентна написанной выше, но выражена гораздо более экстравагантно: In[122]:= Product[i^((-1)^(i-1)), {i,1,2*n}] Out[122]=Gamma[1/2+n]/(Sqrt[Pi]*Gamma[1+n]) Именно в такой форме, конечно, ее и использовал Рамануджан. Однако в случае компьютера мы почему-то твердо уверены, что эта новая фор- мула не внушена богиней Намаккал, а получена в результате честного (и довольно трудного) вычисления — не верю я в эту древнюю восточ- ную мудрость! Если нам трудно представить, что подобное вычисление может быть проведено человеком, то это только потому, что наше соб- ственное математическое воспитание проходило под девизом простоты и неизбежности. И теперь у нас, конечно, не может быть сомнений, что Mathematica су- меет перейти к пределу в предыдущей формуле: In[123]:= Product[i^((-1)^i), {i,1,Infinity}] 116 Out[123]=Sqrt[Pi/2] В действительности, она, конечно, не использует здесь предыдущую фор- мулу, так как вычислить конечное произведение намного сложнее. • Вычисление пределов. Основной командой для вычисления преде- лов в Mathematica является Limit. А именно, для вычисления lim x →c f (x) эта команда вызывается в следующем экзотическом формате: Limit[f[x],x->c]. Обратите внимание, что в угоду традиционному математическому обозна- чению значение параметра c задается так, как будто это опция! В дей- ствительности, конечно, было бы правильно поступать ровно наоборот, а именно, аналистам следовало бы предел функции f точке c записывать не как lim x →c f (x), а просто как lim c (f ), ведь на самом деле здесь никто никуда не переходит, никто ни к кому не стремится, и никакой буковки x в вы- ражении предел функции f точке c нет и отродясь не было. Кроме того, запись f 7→ lim c (f ) подчеркивала бы аналогию с другими гомоморфизмами кольца функций, в которых c выступает в качестве параметра, скажем, с гомоморфизмом эвалюации f 7→ ev c (f ) = f (c). Но чего только не сделаешь ради поддержания самой абсурдной традиции! Однако подобные глупости не мешают системе правильно вычислять пределы, причем в символьном виде: In[124]:= Limit[(Cos[m*x]-Cos[n*x])/(x^2),x->0] Out[124]=1/2*(-m^2+n^2) Команда Limit ищет предел и в тех случаях, когда он не является един- ственным: In[125]:= Limit[Sin[1/x],x->0] Out[125]=Interval[ {-1,1}] Разумеется, к точно такому же результату приведет и вычисление In[126]:= Limit[Sin[x],x->Infinity] Для вычисления левого предела lim x →c− f (x) и правого предела lim x →c+ f (x) ко- манда Limit вызывается в следующих форматах: Limit[f[x],x->c,Direction->1] Limit[f[x],x->c,Direction->-1]. Однако нам встречались случаи, когда система не в состоянии найти предел в символьном виде. В этом случае можно использовать команду NLimit, содержащуюся в пакете NumericalCalculus. 117 § 8. Производные Don’t drink and derive Motto of the society: Mathematicians Against Drunk Deriving В Mathematica имеются две команды символьного дифференцирования D и Derivative, при помощи которых можно явно вычислить все производ- ные, которые могут встретиться любому нематематику. • Дифференцирование. Дифференцирование представляет собой чи- сто алгоритмическую операцию и производная любой элементарной функ- ции может быть вычислена при помощи чисто механических манипуляций. В системе Mathematica имеются две основные команды для вычисления производной, а именно, ◦ функция вычисления производной D, которая находит значение произ- водной; ◦ оператор дифференцирования Derivative[1], в операторной форме ', который вычисляет производную как функцию. В простейшем виде вычисление производной имеет формат D[f[x],x], где f[x] — функция, которую мы хотим продифференцировать, а x — пере- менная, по которой производится дифференцирование. Чтобы убедиться в том, что мы правильно понимаем использование команды D, продифферен- цируем семь функций, производные которых мы хорошо знаем, а именно, функцию x 7→ x n , а также функции exp, ln, cos, sin, arcsin и arccos. Следу- ющее вычисление, как раз, и состоит в нахождении этих семи производных: In[127]:= Map[D[#[x],x]&, {#^n&,Exp,Log,Cos,Sin,ArcSin,ArcCos}] Out[127]= {n*x^(-1+n),E^x,1/x,-Sin[x],Cos[x], 1/Sqrt[1-x^2],-1/Sqrt[1-x^2] } Тому, кто раньше не имел дела с системами компьютерной алгебры, фор- мат этого вычисления может показаться несколько странным, поэтому про- комментируем наиболее важные моменты. Во-первых, мы вычисляем спи- сок производных. Это совершенно необходимо, если мы хотим провести несколько вычислений в одной клетке. Для этого мы приводим список функций, которые хотим продифференцировать, и применяем к элемен- там этого списка дифференцирование D[#[x],x]& при помощи функции Map. Дифференцирование f 7→ f ′ задается здесь в формате анонимной функции, D[#[x],x]&. Напомним, что при этом Slot, в операторной за- писи #, обозначает аргумент чистой функции, в данном случае неизвест- ную функцию от x, которую мы хотим продифференцировать, а Function, в операторной записи &, обозначает применение чистой функции, в данном случае функции дифференцирования по x. Обратите внимание, что в на- шем примере аргумент анонимной функции имеет вид #[x], а не просто #, так что эта неизвестная функция рассматривается именно как функция от того самого x, по которому производится дифференцирование. Может быть чуть проще понять использование этой конструкции во вто- ром примере, когда мы задаем в таком формате функцию возведения в n-ю 118 степень. А именно, с точки зрения Mathematica выражение #^n& является именем функции t 7→ t n в том же самом смысле, в котором sin являет- ся именем функции x 7→ sin(x). Используемый формат имени без явного именования аргументов называется форматом анонимной функции. При некотором навыке использование анонимных функций становится одним из основных средств программирования в Mathematica. Однако начинаю- щему, вероятно, гораздо удобнее использовать формат чистой функции, Function[t,t^n] или, в полной форме Function[t,Power[t,n]], в кото- ром явно указывается, что куда переходит. В этом формате D[#[x],x]& выглядит как Function[f,D[f[x],x]], где явно указано, что функции f сопоставляется ее производная по x. Таким образом, вычисление Function[f,D[f[x],x]][Cos] даст нам - Sin[x], как и просто вычисление D[Cos[x],x]. Чуть более витиеватый формат понадобился нам для того, чтобы вычислить одновременно произ- водные нескольких функций — и то только потому, что мы не знали насто- ящего имени дифференцирования. В действительности, оператор диффе- ренцирования имеет имя, он называется Derivative[1] или, в операторной записи, '. Таким образом, мы могли бы провести то же вычисление с ис- пользованием команды Derivative[1] вместо команды D: In[128]:= Map[#'&, {#^n&,Exp,Log,Cos,Sin,ArcSin,ArcCos}] Out[128]= {n#1^(-1+n)&,E^#1&,1/#1&,-Sin[#1]&, Cos[#1]&,1/Sqrt[1-#1^2]&,-1/Sqrt[1-#1^2]& } Разумеется, первая строчка является просто сокращением строчки In[129]:= Map[Derivative[1], {#^n&,Exp,Log,Cos,Sin,ArcSin,ArcCos}] вычисление которой даст точно такой же результат. Мы видим, что в этом случае вычисляются не значения производных, а сами производные, как функции. Обратите внимание на появляющееся здесь вместо # обозначение аргумента через #1. Дело в том, что если у анонимной функции несколько аргументов, то они обозначаются #1,#2,#3 и так далее. В данном случае программа не знает, в составе какого более сложного выражения мы соби- раемся вызвать эти производные, и нумерует их аргументы с тем, чтобы при появлении следующего аргумента дать ему имя #2 и т.д. • Kettenregel, Produktregel, usw. Разумеется, Mathematica знает все основные правила вычисления производных, причем следующее вы- числение показывает, что она может применять их в символьном виде к неизвестным функциям. Вот как с ее точки зрения выглядят формулы дифференцирования композиции, произведения и частного двух функций, известные под народными названиями Kettenregel, Produktregel и Quo- tientenregel, которые дал им фон Лейбниц: In[130]:= {D[f@g[x],x],D[f[x]*g[x],x],D[f[x]/g[x],x],} Out[130]= {f'[g[x]]*g'[x],g[x]*f'[x]+f[x]*g'[x], f'[x]/g[x]-f[x]*g'[x]/f[x]^2 } Обратите внимание, что f и g обозначают здесь произвольные функции. 119 Разумеется, тем более Mathematica сможет применить эти формулы к кон- кретным функциям. Комментарий. Первая из этих формул — это формула дифференцирования ком- позиции, которая по-русски испокон веку называлась цепным правилом (Ketten- regel, chain rule). К сожалению в последние десятилетия под влиянием украинского языка в русской учебной литературе по так называемой высшей математике укоренился полонизм дифференцирование сложной функции. В польском языке это словообразование вполне оправдано, так как композиция функций называется там сложением (ну а сложе- ние функций, естественно, додаванием). Но в русском языке эта малограмотная калька, впервые появившаяся в украинском переводе учебника Банаха, ничем не подкреплена и ее употребление свидетельствует либо о преступной некомпетентности (если писатель не знает о происхождении этого выражения из польского funkcja zlo ˙zona), либо, в про- тивном случае, о чудовищном манеризме. Кстати, в старых русских учебниках и сама композиция функций называлась по-простому функцией от функции (fonction de fonction), что, конечно, тоже гораздо лучше, чем мифические сложные функции. Следующий диалог иллюстрирует разницу между f −1 и 1/f : In[131]:= {D[InverseFunction[f][x],x],D[1/f[x],x]} Out[131]= {1/f'[f^(-1)[x]],-f'[x]/f[x]^2} При большом желании можно даже выяснить, когда f −1 = 1/f (чрезвы- чайно редко!) Вот, кстати, правило дифференцирования степени (Potenzregel), которое уже далеко не всегда приводится в элементарных учебниках: In[132]:= D[f[x]^g[x],x] Out[132]=f[x]^g[x]*(g[x]*f'[x]/f[x]+Log[f[x]]*g'[x]) Следующий пример показывает, что система действительно полностью отдает себе отчет в том, как применяется цепное правило: In[133]:= {D[ArcCos[Log[x]],x],D[Log[ArcCos[x]],x]} Out[133]= {1/(x*Sqrt[1-Log[x]^2]),-1/(Sqrt[1-x^2]*ArcCos[x])} После этого ясно, что любое другое подобное вычисление является про- сто упражнением в терпении. Вот пример из тех, которые преподаватели высшей математики любят давать на зачетах: In[134]:= D[Log[ArcCos[x]]/ArcCos[Log[x]],x] Out[134]=-1/(Sqrt[1-x&2]*ArcCos[x]*ArcCos[Log[x]])+ Log[ArcCos[x]]/(x*ArcCos[Log[x]]^2*Sqrt[1-Log[x]^2]) • Высшие производные. Понятно, что значение второй производной f ′′ (x) функции f по x можно вычислить при помощи D[D[f[x],x],x], од- нако в языке системы предусмотрены и две более короткие записи этого выражения, а именно D[f[x],x,x] и D[f[x], {x,2}]. В действительности внутренним образом все эти три выражения интерпретируются абсолют- но одинаково, а именно как Derivative[2][f][x], т.е. как значение De- rivative[2][f] в точке x. В свою очередь Derivative[2][f] обозначает собственно вторую производную f ′′ , рассматриваемую как функцию. Ну и, наконец, Derivative[2] обозначает дифференциальный оператор, сопо- ставляющий функции ее вторую производную. Именно для того, чтобы 120 обращаться по имени к этому оператору Mathematica и интерпретирует Derivative как функцию одного аргумента, а именно, порядка!!! У функ- ции Derivative[n] тоже один аргумент, а именно функция, n-ю производ- ную которой мы хотим вычислить. Наконец, у функции Derivative[n][f] снова один аргумент, а именно, точка, в которой мы вычисляем значе- ние f (n) . В Модуле 2 мы подробнейшим образом обсуждаем, чем f[x][y] отличается от f[x,y] и от f[ {x,y}] и почему в разных ситуациях исполь- зуются разные форматы. Однако уже сейчас читатель должен запомнить, что вызов Derivative в формате Derivative[n,f,x] является грубейшей синтаксической ошибкой. Вот, скажем, чему равна вторая производная x x x : In[135]:= Simplify[D[x^x^x,x,x]] Out[135]=x^(-2+x+x^x)(-1+2*x+x^x+x*(3+x+2*x^x)*Log[x]+ x*(2*x+2*x^x+x^(1+x))*Log[x]^2+ x^2*(1+2*x^x)*Log[x]^3+x^(2+x)*Log[x]^4)) Мы применили функцию Simplify, чтобы система привела получающее- ся после применения всех обычных правил выражение для производной к чуть более простому виду, вынеся общий множитель и т.д. Напомним, что система не применяет дистрибутивность автоматически и ей нужно яв- но подсказывать сделать это. Кстати, обратите внимание и на то, что по умолчанию для функции Power используется правая группировка, так что x^y^z истолковывается как x (y z ) , а вовсе не как (x y ) z . Чтобы полу- чить (x y ) z , необходимо явным образом ставить скобки!!! Вычислим для сравнения вторую производную (x x ) x : In[136]:= Simplify[D[(x^x)^x,x,x]] Out[136]=(x^x)^x*(3+2*Log[x]+(x+x*Log[x]+Log[x^x])^2) Теперь уже совершенно ясно, что D[f[x],x,x,x] и D[f[x], {x,3}] выража- ет значение f ′′′ (x) третьей производной f ′′′ функции f по x, и интерпре- тируется как Derivative[3][f][x]. В свою очередь функция f ′′′ вычис- ляется как Derivative[3][f]. Смысл Derivative[n][f] и D[f[x], {x,n}] теперь уже совершенно очевиден. Разумеется, Mathematica умеет применять дифференцирования высших порядков в символьном виде к неизвестным функциям. Вот несколько пер- вых примеров формулы Фаа ди Бруно 37 дифференцирования компози- ции: In[137]:= ColumnForm[NestList[D[#,x]&,f@g[x],4]] Out[137]=f[g[x]] f'[g[x]]*g'[x] g'[x]^2*f''[g[x]]+f'[g[x]]*g''[x] 37 Дональд Кнут предложил называть эту формулу формулой Арбогаста, чем очень огорчил официальный Ватикан. 121 3*g'[x]*f''[g[x]]*g''[x]+g'[x]^3*f'''[g[x]]+ f'[g[x]]*g'''[x] 3*f''[g[x]]*g''[x]^2+6*g'[x]^2*g''[x]*f'''[g[x]]+ 4*g'[x]*f''[g[x]]*g'''[x]+g'[x]^4*f^(4)[g[x]]+ f'[g[x]]*g^(4)[x] Как обычно, f (n) (x) обозначает значение n-й производной d n f dx n • Частные производные. Точно такой же формат как для функций одной переменной имеет и вычисление частных производных по одной из переменных. А именно, D[f[x,y,z], {x,n}] обозначает ∂ n f ∂x n . С другой сто- роны, D[f,x1,x2,x3] обозначает смешанную производную ∂ ∂x 1 ∂ ∂x 2 ∂f ∂x 3 — обратите внимание на порядок переменных!! Напомним, что в такой форме мы обращаемся к значениям этих производных. Собственно сама частная производная ∂ n 1 ∂x n 1 1 ∂ n 2 ∂x n 2 2 ∂ n 3 f ∂x n 3 3 , рассматриваемая как функция, вызывается посредством Derivative[n1,n2,n3][f]. Вот простейший пример In[138]:= {D[Exp[x^2+y^2],x,x],D[Exp[x^2+y^2],x,y], {D[Exp[x^2+y^2],y,y]} Out[138]= {2*E^(x^2+y^2)+4*E^(x^2+y^2)*x^2, 4*E^(x^2+y^2)*x*y, 2*E^(x^2+y^2)+4*E^(x^2+y^2)*y^2 } Как всегда, Mathematica может считать в символьном виде с неизвестными функциями. Вот, скажем, как она понимает цепное правило для функций двух переменных. Следующие формулы можно найти на первых страницах любого учебника по анализу функций нескольких переменных (multivari- able calculus): In[139]:= Simplify[D[f[g[x,y],h[x,y]],x]], Out[139]= {f^(0,1)[g[x,y],h[x,y]]*h^(0,1)[x,y] g^(0,1)[x,y]*f^(1,0)[g[x,y],h[x,y], In[139]:= Simplify[D[f[g[x,y],h[x,y]],y]] Out[139]=f^(1,0)[g[x,y],h[x,y]]*g^(1,0)[x,y] f^(0,1)[g[x,y],h[x,y]]*h^(1,0)[x,y] Как обычно, через f (i,j) (x, y) здесь обозначена функция ∂ i ∂x i ∂ j f ∂x j (x, y). А вот соответствующее правило для вторых производных большинство чита- телей уже вряд ли видело в явном виде. In[140]:= Simplify[D[f[g[x,y],h[x,y]],x,x]] Out[140]=f^(0,2)[g[x,y],h[x,y]]*h^(1,0)[x,y]^2+ 2*g^(1,0)[x,y]*h^(1,0)[x,y]*f^(1,1)[g[x,y],h[x,y]]+ g^(1,0)[x,y]^2*f^(2,0)[g[x,y],h[x,y]]+ f^(1,0)[g[x,y],h[x,y]]*g^(2,0)[x,y]+ f^(0,1)[g[x,y],h[x,y]]*h^(2,0)[x,y] 122 In[141]:= Expand[D[f[g[x,y],h[x,y]],x,y]] Out[141]=h^(0,1)[x,y]*f^(0,2)[g[x,y],h[x,y]]*h^(1,0)[x,y]+ h^(0,1)[x,y]*g^(1,0)[x,y]*f^(1,1)[g[x,y],h[x,y]]+ g^(0,1)[x,y]*h^(1,0)[x,y]*f^(1,1)[g[x,y],h[x,y]]+ f^(1,0)[g[x,y],h[x,y]]*g^(1,1)[x,y]+ f^(0,1)[g[x,y],h[x,y]]*h^(1,1)[x,y]+ g^(0,1)[x,y]*g^(1,0)[x,y]*f^(2,0)[g[x,y],h[x,y]] In[142]:= Simplify[D[f[g[x,y],h[x,y]],y,y]] Out[142]=h^(0,1)[x,y]^2*f^(0,2)[g[x,y],h[x,y]]+ f^(0,1)[g[x,y],h[x,y]]*h^(0,2)[x,y]+ g^(0,2)[x,y]*f^(1,0)[g[x,y],h[x,y]]+ 2*g^(0,1)[x,y]*h^(0,1)[x,y]*f^(1,1)[g[x,y],h[x,y]]+ g^(0,1)[x,y]^2*f^(2,0)[g[x,y],h[x,y]] Явно написать аналог формулы Фаа ди Бруно для высших производных композиции функций нескольких переменнных, ничего при этом не про- пустив, определенно находится за пределами комбинаторных способностей обычного человека. Например, уже в D[f[g[x,y],h[x,y]]x,x,y,y] вхо- дит 46 слагаемых, каждое из которых является произведением двух, трех, четырех или пяти частных производных. • Ряд Тэйлора. Основной командой для вычисления степенных разло- жений функции f является Series. Вызванная в формате Series[f[x], {x,c,m}] эта команда находит первые m членов ряда Тэйлора функции f от x в окрестности точки c. В тех случаях, когда Mathematica может фактически вычислить производные, она вычисляет их, в противном случае — остав- ляет их в символьной форме: In[143]:= Series[f[x], {x,c,4}] Out[143]=f[c]+f'[c]*(x-c)+1/2*f''[c]*(x-c)^2+ 1/6*f'''[c]*(x-c)^3+1/24*f^(4)[c]*(x-c)^4+O[x-c]^5 В тех случаях, когда команда Series не может найти ряд Тэйлора (напри- мер, если значения каких-то производных в точке c бесконечны), она пишет разложение функции в ряд, в который входят также отрицательные и/или дробные степени x − с: In[144]:= Series[Exp[Sqrt[x]], {x,0,3}] Out[144]=1+x^(1/2)+1/2*x+1/6*x^(3/2)+1/24*x^2+ 1/120*x^(5/2)+1/720*x^3+ O[x]^(7/2) 123 § 9. Интегралы Are you an Analyst? — I am an Oralist. George Bergmann 3 √ 3 ∫ 1 z 2 dz cos(3π/9) = ln( 3 √ 3) Anonimous Analyst 38 Одной из самых мощных команд системы является команда Integrate, при помощи которой вычисляются как неопределенные, так и определен- ные интегралы от функций одной и нескольких переменных. Обратим вни- мание, что имплементация команды Integrate включает 600 страниц кода на C и еще около 500 страниц высокоуровневого кода! • Неопределенный интеграл. Первообразная функции f вычисляет- ся при помощи команды Integrate, вызываемой в том же формате, что и команда дифференцирования D, а именно, Integrate[f[x],x]. Эта ко- манда возвращает какую-то первообразную функции f . Таким образом, взятие неопределенного интеграла считается операцией обратной диффе- ренцированию: In[145]:= {D[Integrate[f[x],x],x],Integrate[D[f[x],x],x]} Out[145]= {f[x],f[x]} Вот пара простеньких интегралов, которые легко берутся интегрированием по частям: In[146]:= {Integrate[ArcTan[x],x],Integrate[ArcSin[x],x]} Out[146]= {x*ArcTan[x]-1/2*Log[1+x^2],Sqrt[1-x^2]+x*ArcSin[x]} Конечно, система может легко взять любой стандартный интеграл, для вы- числения которого не требуется ничего, кроме многократного применения обычных формул для элементарных функций и стандартных подстановок: In[147]:= Simplify[Integrate[Cos[x]^10,x]] Out[147]=1/10240*(2520*x+2100*Sin[2*x]+600*Sin[4*x]+ 150*Sin[6*x]+25*Sin[8*x]+2*Sin[10*x]) In[148]:= Integrate[(x+Sqrt[x^2+1])/(x+1-Sqrt[x^2+1]),x] Out[148]=1/2*(x+x^2+Sqrt[1+x^2]+x*Sqrt[1+x^2]+ ArcSinh[x]+2*Log[x]-Log[1+Sqrt[1+x^2]) Появление гиперболического арксинуса во втором интеграле совершенно естественно, если вспомнить, что Z 1 √ 1 + x 2 dx = arcsh(x) + c. Тем не менее, выводить вручную даже такие совсем простые формулы не слишком приятно. 38 Для удобства оралистов приведем американское чтение этого лимерика: Integral z-squared dz/ from 1 to the cube root of 3/ times the cosine/ of three π over 9/ equals log of the cube root of e. 124 • Неберущиеся интегралы. В отличие от вычисления производных вычисление первообразных представляет собой творческое занятие. Тон- кий момент здесь состоит в том, что интеграл от элементарной функции уже совершенно не обязательно является элементарной функцией, но ни- каких очевидных критериев, которые позволяют сказать, когда подобная неприятность случается, нет. Уже композиция двух основных элементар- ных функций может не иметь элементарного интеграла. Вот несколько совсем простеньких примеров, которые хорошо иллюстрируют возникаю- щие здесь проблемы. ◦ Интеграл R sin(ln(x)) dx легко берется интегрированием по частям. Каждый из нас должен был проделать следующее вычисление на первом курсе: In[149]:= Integrate[Sin[Log[x]],x] Out[149]=(1/2)*x*Cos[Log[x]]+(1/2)*x*Sin[Log[x]] ◦ Но вот переставив здесь Sin и Log, мы получим интеграл R ln(sin(x)) dx, в который входит полилогарифм: In[150]:= Integrate[Log[Sin[x]],x] Out[150]=(-x)*Log[1-E^(2*I*x)]+x*Log[Sin[x]]+ (1/2)*I(x^2+PolyLog[2,E^(2*I*x)]) где PolyLog[n,x] обозначает полилогарифм Li n (x), который определяет- ся как Li n (x) = ∞ X k=1 z k k n . ◦ Также и функция ln 2 (x) = ln(ln(x)) не интегрируется в элементарных функциях. А именно, интеграл R ln 2 (x) dx равен In[151]:= Integrate[Log[Log[x]],x] Out[151]=x*Log[Log[x]]-LogIntegral[x] Здесь LogIntegral[x] обозначает интегральный логарифм li(x), кото- рый можно определить как li(x) = Z x 0 dt ln(t) . ◦ А вот R sin 2 (x) dx вообще никак не выражается в терминах более из- вестных функций, так что если Вы предложите Mathematica вычислить In[152]:= Integrate[Sin[Sin[x]],x] она после некоторого раздумья просто перепишет это выражение. Это зна- чит, что системе не известно никакой более простой или более явной формы для Integrate[Sin[Sin[x]],x]. Мы не сомневаемся, что к этому моменту читатель уже в состоянии отличить неберущийся интеграл R sin 2 (x) dx от интеграла R sin(x) 2 dx = x/2 − sin(2x)/4. 125 ◦ Чтобы завершить тему с возведением в квадрат всего, что в sin(x) возводится в квадрат, вычислим еще интеграл R sin(x 2 ) dx. С точностью до нормировки ответ называется синус интегралом Френеля: In[153]:= Integrate[Sin[x^2],x] Out[153]=Sqrt[Pi/2]*FresnelS[Sqrt[2/Pi]*x] Этот интеграл, как и косинус интеграл Френеля R cos(x 2 ) dx или, после нормировки, R cos(πx 2 /2) dx, часто встречается в оптике. Да чего там композиции — уже произведения двух основных элементар- ных функций совершенно не обязательно интегрируются в элементарных функциях (произведение двух производных вообще не обязано быть произ- водной чего бы то ни было). Вот два совсем простых примера. При попытке вычислить Z x sin(x) dx встречается уже знакомый нам полилогарифм In[154]:= Integrate[x/Sin[x],x] Out[154]=x*(Log[1-E^(I*x)]-Log[1+E^(I*x)])+ I*(PolyLog[2,-E^(I*x)]-PolyLog[2,E^(I*x)]) А вот Z sin(x) x dx имеет специальное название – интегральный синус: In[155]:= Integrate[Sin[x]/x,x] Out[155]=SinIntegral[x] Через SinIntegral[x] и CosIntegral[x] в Mathematica как раз и обозна- чаются интегральный синус и интегральный косинус соответственно: Si(x) = Z x 0 sin(t) t dt, Сi(x) = − Z ∞ x cos(t) t dt, Каждый может узнать все эти факты секунд за 30, со скоростью нуж- ной, чтобы набрать эти интегралы на клавиатуре. Заметим, что большая часть этих интегралов вообще не упоминается ни в общем курсе матема- тического анализа, ни в большинстве стандартных учебников. Впрочем, и в этом отношении замечательное руководство Владимира Антоновича Зо- рича 39 оказывается с огромным отрывом лучшим из всех учебных текстов на русском языке. А именно, в упражнениях 6 и 7 на стр. 380–381 перво- го тома не только перечисляются все эти интегралы, но и объясняется их происхождение из интегральной экспоненты: Ei(x) = Z x −∞ e t t dt, на языке Mathematica, естественно, ExpIntegralEi[x]. Поэтому, если Вы не являетесь специалистом по некоторым разделам анализа или математической физики, и не читали Зорича, то Вы, скорее 39 В.А.Зорич, Математический анализ. т. I,II. — М., МЦНМО, 2002, с.1–657, с.1–787. 126 всего, вообще не видели большинства этих интегралов. А теперь скажите, сколько времени Вам понадобилось бы, чтобы извлечь эту информацию из традиционного справочника, ну скажем, из 40 ? • Определенный интеграл. Определенный интеграл R b a f (x) dx вы- числяется при помощи той же самой команды Integrate, которая в этом случае вызывается в формате Integrate[f[x], {x,a,b}], где явно указа- ны пределы интегрирования. Mathematica знает, как дифференцировать определенный интеграл по пределам интегрирования: In[156]:= D[Integrate[f[t], {t,g[x],h[x]}],x] Out[156]=-f[g[x]]*g'[x]+f[h[x]]*h'[x] А вот формула Ньютона—Лейбница: In[157]:= f[x ]:=D[g[x],x]; Integrate[f[x], {x,a,b}] Out[157]=-g[a]+g[b] Еще раз объясним, что здесь происходит. В строке ввода две команды, раз- деленные точкой с запятой. Первая из них — это операция отложенного присваивания := при помощи которой мы задаем равенство f = g ′ . Как всегда при этом бланк используется для того, чтобы сообщить системе, что f (x) = g ′ (x) для всех x. Вместо этого мы могли при желании про- извести немедленное присваивание f=Derivative[1][g], но в этом случае нужно задавать саму функцию f , а не ее значения! После этого вторая команда предлагает системе проинтегрировать f . Только что сказанное убеждает нас в том, что система сможет взять любой определенный интеграл, если она может взять соответствующий неопределенный интеграл и значения первообразной на концах конечны. Гораздо интереснее, как она будет выкручиваться при возникновении неопределенностей или бесконечных значений! Оказывается, это тоже не представляет проблемы: In[158]:= Integrate[Log[x], {x,0,1}] Out[158]=-1 Мы уже знаем, что интеграл от функции sin 2 (x) = sin(sin(x)) не берется. Тем не менее, некоторые определенные интегралы от этой функции берутся явно: In[159]:= Integrate[Sin[Sin[x]], {x,0,Pi}] Out[159]=Pi*StruveH[0,1] Появляющаяся в этой формуле функция Струве StruveH является спе- циальным решением неоднородного уравнения Бесселя. Вот еще один не самый банальный интеграл, вычисление которого обыч- но сводится к многомерным интегралам, либо интегралам, зависящим от параметра: 40 А.П.Прудников, Ю.А.Брычков, О.И.Маричев, Интегралы и ряды. Элементарные функции. — Наука, 1981. с.1—798. 127 In[160]:= Integrate[ArcTan[x]/(x*Sqrt[1-x^2]), {x,0,1}] Out[160]=1/2*Pi*ArcSinh[1] По форме ответа кажется, что Mathematica в лоб вычислила соответ- ствующий неопределенный интеграл и использовала формулу Ньютона— Лейбница. Однако это абсолютно не так, поскольку вычислять неопреде- ленный интеграл Z arctg(x) x √ 1 − x 2 dx она не умеет 41 !!! Таким образом, для нас остается загадкой, каким образом Mathematica получила правильный от- вет в данном случае! Впрочем, даже зная этот ответ, неспециалист, скорее всего, далеко не сразу заметит, что это в точности π 2 ln(1 + √ 2). • Несобственные интегралы. Одними из самых знаменитых и часто встречающихся в приложениях определенных интегралов являются инте- грал Эйлера—Пуассона, известный также как Гауссов интеграл 42 , и интеграл Дирихле: Z ∞ −∞ e −x 2 dx = √ π, Z ∞ −∞ sin(x) x dx = π. Посмотрим, может ли Mathematica взять эти интегралы, ведь пределы ин- тегрирования в них бесконечны. Оказывается, и это не представляет ника- кой проблемы: In[161]:= Integrate[E^(-x^2), {x,-Infinity,Infinity}] Out[161]=Sqrt[Pi] In[162]:= Integrate[Sin[x]/x, {x,-Infinity,Infinity}] Out[162]=Pi Вот еще один очень красивый определенный интеграл, соответствующий неопределенный интеграл выражается в терминах полилогарифма, но при интегрировании от 1 до ∞ все исчезает, остается (пространство, звезды и певец) константа Каталана C: In[163]:= Integrate[Log[x]/(1+x^2), {x,1,Infinity}] Out[163]=Catalan Заметим, что Z 1 0 ln(x) 1 + x 2 dx = −C. Если хотите испытать культурный шок, то попробуйте теперь вычислить Z ∞ 0 ln(x) 1 + x 2 dx. 41 Это один из немногих известных нам интегралов, которые Mathematica не умеет вычислять. Впрочем, в подобных случаях больной все равно должен обращаться к док- тору, поскольку на стр. 268 книги “Интегралы и ряды” этого интеграла тоже нет, а есть только более простой интеграл ∫ x arctg(x) √ 1 − x 2 dx. Упражнение для профессионалов: зная ответ придумайте подстановку, вычисляющую этот интеграл. 42 В теории вероятностей интеграл erf(x) = ∫ x 0 e −t 2 dt, рассматриваемый как функция от верхнего предела, называется еще функцией ошибок или интегралом ошибок = error function. 128 Попытка вычислить интеграл Integrate[1/x, {x,-1,1}] не приведет к большому успеху, так как система выдаст следующее сообщение об ошибке Integrate: Integral of 1/x does not converge on {-1,1}. В подобном случае можно попытаться вычислить главное значение инте- грала. Вычисление In[164]:= Integrate[1/x, {x,-1,1},PrincipalValue->True] даст значение 0. Правило PrincipalValue->True меняет установку одной из опций функции Integrate. Это типичный пример того, как настройка опций меняет способ работы функций и влияет на получающийся ответ. • Кратные интегралы. Команда Integrate служит и для вычисления кратных интегралов. Скажем, при вычислении двойного интеграла она вызывается в формате Integrate[f[x,y], {x,a,b},{y,c,d}]. Обратите внимание, что, как и в командах наподобие Table или Array, порядок итераторов здесь таков: внутренние итераторы изображаются по- следними. Таким образом, Integrate[f[x,y], {x,a,b},{y,c,d}] это Z b a Z d c f (x, y) dy dx, а вовсе не Z d c Z b a f (x, y) dx dy. Заметим, что при вычислении кратных интегралов результат в большин- стве случаев следует упрощать, так как без этого система возвращает его as is, как сумму большого количества слагаемых, в том виде, как они возникают. Вот один из классических двойных интегралов, In[165]:= Simplify[Integrate[1/Sqrt[x^2+y^2], {x,-1,1},{y,-1,1}]] Out[165]=2*(2*ArcSinh[1]-Log[-1+Sqrt[2]]+Log[1+Sqrt[2]]) В этот момент у каждого, кто когда-либо преподавал Multivariable Cal- culus, да и просто у того, кто с самого начала внимательно читал этот параграф, закрадываются подозрения, что это все еще не окончательный ответ. Еще одна попытка заканчивается полной победой: In[166]:= FullSimplify[Integrate[1/Sqrt[x^2+y^2], {x,-1,1},{y,-1,1}]] Out[166]=8*ArcSinh[1] Математик — это тот, кто никогда не останавливается на достигнутом, од- нако такого тройного интеграла, пожалуй, не считали даже наши студенты: In[167]:= FullSimplify[Integrate[1/Sqrt[x^2+y^2+z^2], {x,-1,1},{y,-1,1},{z,-1,1}]] Out[167]=-2*(Pi-4*ArcCsch[Sqrt[2]]+Log[-8(-2+Sqrt[3])]- 6*Log[1+Sqrt[3]]) Естественно, пределы интегрирования сами могут быть функциями от дру- гих аргументов: In[168]:= Integrate[1, {x,-1,1},{y,-Sqrt[1-x^2],Sqrt[1-x^2]}, {z,-Sqrt[1-x^2-y^2],Sqrt[1-x^2-y^2]}] |