Математические эксперименты с Mathematica

Преподаватель : Тиморин Владлен Анатольевич
Кафедра геометрии и топологии (комната 314)
Адрес электронной почты : vtimorin@hotmail.com
Присутственные часы : четверг 14-15:30, пятница 15 : 30 - 16 : 50

Внизу находится список упражнений, которые нужно сдать на занятиях.

Учебные материалы, которые здесь имеются, находятся в стадии разработки - следите за обновлениями.

1. Введение

Цель курса

Пользуясь системой Mathematica, мы будем ставить математические эксперименты. Попутно (это не является главной целью) мы освоим основные принципы работы с системами компьютерной алгебры и основные принципы программирования на языках высокого уровня. В курсе будет много примеров из  различных разделов математики - на этих примерах можно увидеть и почувствовать различные математические объекты (такие, как функции, перестановки, множества, графы и т.д.).

Система оценок

Итоговая оценка складывается из оценок за работу в классе и оценок за индивидуальные проекты.

Что такое Mathematica?

Mathematica - это одновременно система компьютерной алгебры и язык программирования. Система компьютерной алгебры - это, грубо говоря, очень продвинутый калькулятор, позволяющий производить не только и не столько численные, сколько символические вычисления - раскладывать многочлены на множители, считать производные и интегралы, работать с векторами и матрицами, и т.д.
Как язык программирования, система Mathematica включает в себя все основные возможности таких языков, как C++. Однако, в отличие, например, от C++, многие стандартные операции (скажем, динамическое выделение памяти) происходят "за кадром". Поэтому типичная программа в системе Mathematica гораздо короче и прозрачней, чем аналогичные программы на языках программирования более низкого уровня.

Разработкой системы Mathematica занимается компания Wolfram Research. Основным коммерческим конкурентом системы Mathematica является программа Maple, созданная канадским университетом Ватерлоо. Есть также некоммерческие (в основном, бесплатные и свободно распространяемые) системы компьютерной алгебры. К ним относится, например, Yacas (аббр. от yet another computer algebra system). Сейчас большинство работающих математиков, физиков, инженеров пользуются той или иной системой компьютерной алгебры. Какую систему выбрать - вопрос вкуса и доступности (например, не всякому по карману приобретение индивидуальной лицензии системы Mathematica). В любом случае, предлагаемый курс облегчит освоение не только системы Mathematica, но и любой другой аналогичной программы.

Как пользоваться системой Mathematica?

Система Mathematica состоит из ядра (т.е. программы, производящей вычисления) и оболочки (которая используется для ввода запросов, вывода результатов, а также для текстовых комментариев).

Для того, чтобы система Mathematica произвела вычисление, нужно ввести задание (последовательность команд) и нажать одновременно клавиши Shift и Enter.

Оболочка системы Mathematica является, по существу, полноценным текстовым редактором. В ней можно выбирать шрифты, форматы, набирать специальные символы и математические формулы. Есть даже книги, написанные целиком в системе Mathematica.

2. Функции и графики

Арифметические операции

Конечно, можно использовать систему Mathematica как калькулятор.

2+2

2 3

2^10

100!

2/(3+2/3)

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

Sqrt[12]

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

2.^(1/2)

N[Sqrt[2],100]

N[%^2,200]

На карманном калькуляторе этого не видно! Число (e^(π(Sqrt[163])))целое с точностью до 12 знаков после запятой. Число 163 замечательно тем, что это самое большое число d,
для которого кольцо целых поля Q[(Sqrt[-d])] является кольцом главных идеалов.

N[E^3,50]

Элементарные функции

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

Sin[Pi/6-.00001]

Sin[Pi/5]

N[ArcSin[Pi]]

N[Sin[Pi/7]]

N[Exp[2+3 I],5]

ArcTan[1]

Составим шпаргалку специальных значений синуса (пригодится абитуриентам).

X={0,1/6,1/5,1/3,2/5,1/2};
Y=Map[Sin[# Pi]&,X];
Grid[{Pi X,Y},Frame->All,Background->LightBlue]

Специальные функции

Гамма - функция

Замечательное наблюдение Эйлера состоит в том, что функция Γ(x)=(x-1)!, определенная комбинаторно для целых положительных значений x, является ограничением некоторой аналитической функции.

Gamma[6]

Gamma[1/2]

Вот график гамма - функции для вещественных значений x.

Plot[Gamma[x],{x,-5,5}]

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

С другой стороны, гамма - функция определена и для комплексных значений параметра. Чтобы увидеть, где у нее особенности на комплексной плоскости, нарисуем график модуля.

Plot3D[Abs[Gamma[x+I y]],{x,-5,5},{y,-2,2},MaxRecursion->5]

Отсюда видно, что никаких других особенностей, кроме уже указанных, нет.

Основное свойство гамма - функции, обобщающее соответствующее свойство факториала, таково : Γ(s+1)=sΓ(s).

FullSimplify[s Gamma[s]]

Интересны явные формулы для гамма - функции.

(Gamma[s])=(lim_(n->∞)) ((n! (n^s))/((∏_(k=0)^n)(k+s)))

Limit[(n!*n^s)/Product[k + s, {k, 0, n}], 
   n -> Infinity];
FullSimplify[%]

(Gamma[s])=(∫_0^∞)exp(-t) (t^(s-1))dt            если    Re(s)>0.

Integrate[Exp[-t]*t^(s-1),{t,0,Infinity},Assumptions->Re[s]>0]

Функциональное уравнение для гамма - функции :

(Gamma[s]) (Gamma[1-s])=(π/(sin(π s)))

FullSimplify[Gamma[s] Gamma[1-s]]

Дзета - функция

Дзета-функция Римана очень важна в теории чисел. Она определяется так:

(Zeta[s])=(∑_(n=1)^∞)(1/(n^s))

Вот некоторые специальные значения дзета - функции :

Manipulate[Zeta[k],{k,-5,5,1},{k,Range[-5,5]}]

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

(Zeta[s])=(∏_(n=1)^∞)(1/(1-((((Prime[n])))^(-s))))

Product[(1-Prime[n]^(-s))^(-1),{n,1,Infinity}]

Интересно расположение нулей дзета - функции.

ContourPlot[Abs[1/Zeta[x+I y]],{x,-2,2},{y,-40,40},Contours->{0.5,1,2,5,10^10},AspectRatio->2,MaxRecursion->3]

На этом графике можно лучше увидеть нули дзета - функции.

Plot[Abs[Zeta[1/2+I x]],{x,-40,40},Filling->Axis,ColorFunction->Function[{x,y},Hue[y]]]

Функциональное уравнение для дзета-функции:

(((Zeta[s]) (Gamma[(s/2)]))/(π^(s/2)))=(((Zeta[1-s]) (Gamma[((1-s)/2)]))/(π^((1-s)/2)))

Полюс в 1

Limit[Zeta[s]-1/(s-1),s->1]

N[%]

Переменные и функции

x=5

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

y=5;z=5;

{x,y}

x==y

x+2

Clear[x]
x+2

Clear[f]
f[x_]:=x^2

f[2]

D[f[x],x]

Clear[g]
g:=#^2&
g[2]

Можно давать рекуррентные определения :

Clear[fib];
fib[0]:=1; fib[1]:=1;
fib[i_]:=fib[i]=fib[i-1]+fib[i-2]

fib[5]

Манипулятор

Manipulate[Plot[Sin[n x],{x,0,6}],{n,1,10}]

Manipulate[Factor[x^n-1],{n,10,100,1}]

Manipulate[Plot[Sin[a x+b],{x,0,6}],{{a,2,"Multiplier"},1,4},{{b,0,"Phase Parameter"},0,10}]

Manipulate[Nest[Subsuperscript[#,#,#]&,x,n],{n,1,6,1}]

Manipulate[{x^2,Button["reset",x=0]},{x,0,10}]

Manipulate[Plot[f[x],{x,0,2Pi}],{f,{Sin,Cos,Tan,Cot}}]

Manipulate[Graphics[{color,Disk[]}],{color,Purple}]

Manipulate[Graphics[Polygon[pt],PlotRange->2],{{pt,{{0,0},{1,0},{1,1},{0,1},{1,-1}}},Locator}]

3. Алгебра

Упрощение

Simplify[x^2+2 x y +y^2]

Simplify[Sqrt[x^2]]

Simplify[Sqrt[x^2],Assumptions->x>0]

Simplify[Log[Sin[x]^2+Cos[x]^2]]

Операции над многочленами

Коэффициент многочлена (((x+1))^3) при второй степени x:

Coefficient[(x+1)^3,x,2]

Coefficient[(x+y)^4,x y^3]

CoefficientList[(x+3)^5,x, Modulus->7]

Expand[(x+2)^3 (x-1)^2]

Simplify[(x-1)(x+1)(x^2+1)+1]

PolynomialGCD[(1+x)^2(2+x)(4+x),(1+x)(2+x)(3+x)]

Производящие функции

Производящая функция для числа разбиений (=диаграмм Юнга)

F:=Product[1/(1-t^k),{k,1,Infinity}]
plist[n_]:=CoefficientList[Series[F,{t,0,n}],t]
plist[20]

ListPlot[plist[100]//Log]

Задача 3.1. Экспериментально, найдите асимптотику для p (n) = количества разбиений натурального числа n.

Числа Бернулли

Clear[t];
F=t/(Exp[t]-1);
b[n_]:=Coefficient[Series[F,{t,0,50}],t,n] n!;
Table[b[n],{n,1,10}]
Table[BernoulliB[n],{n,1,10}]

Задача 3.2. Исследуйте асимптотику чисел Бернулли с четными номерами.

Многочлены Бернулли и формула Эйлера-Маклорена

Формула Эйлера - Маклорена выражает связь между суммой ряда

(∑_(n=1)^∞)f(x+n)

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

F(x)=(∫_x^∞)f(t)dt

А именно :

(∑_(n=1)^∞)f(n+x)=((D )/((e^D)-1))F(x),

где D обозначает дифференцирование по x.

Например (хотя, очевидно, формула Эйлера - Маклорена не применима к данному случаю):

(∑_(n=1)^∞)(((n+x))^k)=((D )/((e^D)-1))(((x^(k+1))/(k+1)))

Многочлен в правой части называется многочленом Бернулли (Subscript[B, k])(x).

Table[BernoulliB[n,z],{n,0,5}]

Задача 3.3. Выразите сумму

(∑_(i=1)^n)(i^k)

через многочлены Бернулли.

Числа Бернулли (Subscript[B, k]) - это значения свободных коэффициентов в многочленах Бернулли.

Table[BernoulliB[n,0],{n,0,10}]
Table[BernoulliB[n],{n,0,10}]

Задача 3.4. В многочлене (((B+x))^n), замените все степени переменной B соответствующими числами Берлулли, т.е. (B^k) заменяется на (Subscript[B, k]). Проверьте, что получится многочлен Бернулли (Subscript[B, n])(x). Проделайте это вычисление для некоторых n, скажем, для n=18,19,20.

Задача 3.5. Рассмотрим функцию

((t exp(t x))/(exp(t)-1))

Найдите коэффициенты разложения этой функции по степеням t (скажем, первые 20 коэффициентов). Проверьте, что коэффициенты совпадают с многочленами Бернулли.

Задача 3.6. На основании экспрериментальных данных, угадайте формулы для явных значений дзета - функции ζ (-n) и ζ (2 n).

Разложение на множители, элементарные дроби

Factor[8-4 x-10 (x^2)+(x^3)+4 (x^4)+(x^5)]

Factor[x^10-1,Modulus->2]

Factor[x^2-2]

Factor[x^2-2, Extension -> Sqrt[2]]

Factor[1+x^4,Extension->{Sqrt[2],I}]

Apart[1/((1+x)(5+x))]

Apart[1/((x+a) (x+b) (x+c)),x]

Дискриминанты и результанты

Рассмотрим многочлен от 1 переменной, старший коэффициент которого равен 1. Дискриминант такого многочлена определяется как произведение квадратов всех попарных разностей корней.

Discriminant[Product[x-(Subscript[t, j]),{j,1,5}],x]//Simplify

Задача 3.7. Посчитайте дискриминант многочлена (x^n)+1 как функцию от n (на основании экспериментов, угадайте явную формулу). Попробуйте доказать эту формулу.

Задача 3.8. Как определить дискриминант в общем случае (т.е. когда старший коэффициент не обязательно равен 1)?

Discriminant[a x^2+b x+c,x]

Discriminant[x^10-5x^7-3 x+9,x]

Discriminant[a x^3+b x^2+c x +d,x]

Задача 3.9. Особенность "ласточкин хвост" возникает на поверхности, заданной в трехмерном пространстве с координатами a,b,c уравнением
Discriminant[a (x^3)+b (x^2)+c x+(x^5),x]==0
Нарисуйте эту поверхность так, чтобы особенность в начале координат была хорошо видна.

Результант двух многочленов P и Q - это произведение (u-v) по всем парам (u,v), таким, что P(u)=Q(v)=0.

Resultant[(x-a)(x-b),(x-c)(x-d)(x-e),x]

P:=x^2+a x+b
Q:=x^3+c x^2+d x+e
Resultant[P,Q,x]

Результант можно получить, исключая переменную x из системы уравнений P(x)=Q(x)=0.

el=Eliminate[{P==0,Q==0},x]
Simplify[el[[1]]-el[[2]]-Resultant[P,Q,x]]

Задача 3.11. Найдите связь между дискриминантом многочлена и результантом многочлена и его производной (на примере многочленов степени 5).

Результант можно также посчитать как определитель матрицы Сильвестра.

SylvesterMatrix1[poly1_,poly2_,var_]:=Function[{coeffs1,coeffs2},With[{l1=Length[coeffs1],l2=Length[coeffs2]},Join[NestList[RotateRight,PadRight[coeffs1,l1+l2-2],l2-2],NestList[RotateRight,PadRight[coeffs2,l1+l2-2],l1-2]]]][Reverse[CoefficientList[poly1,var]],Reverse[CoefficientList[poly2,var]]]

S:=SylvesterMatrix1[P,Q,x]
MatrixForm[S]

Det[S]-Resultant[P,Q,x]

Решение алгебраических уравнений

ClearAll;

Solve[x^2+a x+1==0,x]

Solve[{a x+y==7,b x-y==1},{x,y}]

Solve[x^3+a x+b==0,x]

Solve[x^5==1,x]

NSolve[x^5==1,x]

Решение уравнения (z^7)=1 через кубические и квадратные радикалы.

z=(-1)^(1/7);
p=MinimalPolynomial[z+1/z];
Solve[p[u]==0,u]

Задача. Какие правильные n-угольники можно построить циркулем и линейкой? Скажем, для n от 1 до 20.

3 1/2. Линейная алгебра

Операции над матрицами

A:={{1,2,3},{-2,1,-1},{1,0,1}}
MatrixForm[A]

MatrixForm[Transpose[A]]

B:=Array[1+#1^2+#2^2&,{5,5}]
MatrixForm[B]

X:=Table[Subscript[x,i],{i,5}]
MatrixForm[X]

MatrixForm[A+B]
MatrixForm[A.B]

Det[A+B]
Tr[B]

Общее выражение для определителя матрицы 3 на 3.

Det[Array[(Subscript[a, #1,#2])&,{3,3}]]

Определитель Вандермонда

MatrixForm[Array[((Subscript[t, #1])^(#2-1))&,{4,4}]]
Det[Array[((Subscript[t, #1])^(#2-1))&,{4,4}]]//Factor

Нормальные формы матриц

Можно приводить матрицы к жордановой нормальной форме.

A={{27,48,81},{-6,0,0},{1,0,3}};
A//MatrixForm
JordanDecomposition[A][[2]]//MatrixForm

Собственные векторы и собственные значения

ListPlot[Eigenvalues[Table[N[1/(i+j+1)],{i,50},{j,50}]]//Log]

Приведение квадратичной формы к диагональному виду.

q = 20 x^2 - 16 x y + 23 y^2 - 12 x z - 2 y z + 17 z^2;
m=Normal[CoefficientArrays[q,{x,y,z},Symmetric->True][[3]]];
vals=N[Eigenvalues[m]]

Цепи Маркова и случайные блуждания

Пусть X - конечное множество состояний (для простоты будем считать, что X - это натуральные числа от 1 до n). Допустим, что последовательность случайных величин принимает значения во множестве X, и вероятность того, что следующая
случайная величина принимает значение i при условии, что предыдущая случайная величина принимает значение j, равна (Subscript[p, ij]).

Матрица P, составленная из вероятностей (Subscript[p, ij]), называется матрицей перехода. Заметим, что вектор распределения вероятностей на каждом следующем шаге умножается на матрицу P.

Инвариантное распределение вероятностей (инвариантная мера на X), если существует, соответствует собственному вектору матрицы перехода с собственным значением 1. Можно показать, что все координаты такого вектора будут неотрицательны. Вектор нужно нормировать таким образом, чтобы сумма координат была равна 1.

invm[A_]:=Module[{P,vals,vecs,i,v},
P=Transpose[A];
{vals,vecs}=Eigensystem[P];
For[i=1,i≤Length[vals],i++,
If[vals[[i]]==1.,
v = vecs[[i]];
Return[v/Sum[v[[k]],{k,1,Length[v]}]]
]
]
]

Случайное блуждание на графе

GraphPlot[{1->2,2->3,3->4,4->5}]

invm[{{0,1.,0,0,0},{.5,0,.5,0,0},{0,.5,0,.5,0},{0,0,.5,0,.5},{0,0,0,1.,0}}]

4. Анализ

Дифференцирование

ClearAll[x,f,g];

Команда D[expr,var] дифференцирует выражение expr относительно переменной var (все остальные входящие в это выражение переменные предполагаются не зависящими от var).

D[x^n,x]

Команда D[expr, var,n] дифференцирует выражение expr относительно переменной var n раз.

D[Sin[x]^10,{x,4}]

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

D[1/(x^2+y^2),x,y]

Производные можно брать в символическом виде. Например, вот как можно продифференцировать произведение двух функций 3 раза по правилу Лейбница:

D[f[x] g[x],{x,3}]

А это многократное применение правила дифференцирования сложной функции и правила Лейбница для произведения функций:

D[f[g[x]],{x,4}]

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

D[x^2+y^2,{{x,y}}]

Напомним, что кривизна кривой x=x(t) (где x - это вектор-функция) может быть вычислена по формуле ((|v×a|)/(|v((|)^3))), где v, a - векторы скорости и ускорения, соответственно.
Кривизна круговой спирали:

Clear[r,x,v,a,b,c]
$Assumptions= r>0 $$ t>0 $$ b>0 $$ c>0;
x={r Cos[b t],r Sin[b t],c t};
v = D[x,t];
a=D[v,t];
FullSimplify[Norm[Cross[v,a]]/Norm[v]^3]

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

Clear[x,a,b]

Можно находить неопределенные интегралы

Integrate[Log[x],x]

Mathematica может явно находить значения многих определенных интегралов, в том числе несобственных:

Integrate[Exp[-x^2],{x,0,Infinity}]

Integrate[1/Sqrt[1+x^2],{x,0,1}]

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

Integrate[1/Sqrt[1+x^2],{x,a,b}]

Integrate[1/Sqrt[1+x^2],{x,a,b},Assumptions->{a>0,b>a}]

Можно сразу вычислять численное значение определенного интеграла.

NIntegrate[1/Sqrt[1+x^4],{x,0,1}]

Чтобы посчитать число нулей функции f в области, ограниченной некоторой кривой на комплексной плоскости, нужно проинтегрировать функцию log(f) вдоль этой кривой, и разделить на 2πi.

Clear[x];
NIntegrate[(D[Log[Sin[x]],x]x)/.x->30 Exp[2 Pi I t],{t,0,1}]

Числовые ряды

Mathematica знает, как вычислять многие специальные суммы и произведения:

Sum[i^2,{i,1,n}]

Это относится и к бесконечным суммам, т.е. к рядам :

Sum[1/n^2,{n,1,Infinity}]

При вычислении сумм рядов, интегралов и проч. используются значения специальных функций. Например:

Sum[1/n^3,{n,1,Infinity}]

Mathematica распознает многие замечательные константы, например, константу Эйлера-Маскерони:

Sum[1/n-Log[(n+1)/n],{n,1,Infinity}]

Sum[1/(1+i+j)^3,{i,1,Infinity},{j,1,i}]

Ряды Тейлора, Фурье, Пьюзо

f=Series[Exp[x],{x,0,10}]

Normal[f]

g=Series[Exp[-x],{x,0,10}]
Log[g]

f g
f/.x->g-1

h:=Series[Log[1+x],{x,0,10}]
h/.x->f-1
ComposeSeries[h,f-1]

InverseSeries[f-1]

Разложение неявных функций в ряды Пьюзо.

ClearAll;
Series[y//.Solve[x^3+x y+y^4==0,y],{x,0,5}]

Ортогональные многочлены

Table[LegendreP[n,x],{n,1,10}]

p[n_]:=Simplify[1/(2^n n!) D[(x^2-1)^n,{x,n}]]
Table[p[n],{n,1,10}]

Integrate[p[8] p[9],{x,-1,1}]

de[n_,f_]:=(1-x^2)D[f,{x,2}]-2x D[f,x]+n(n+1)f

de[8,p[8]]//Simplify

Plot[Table[LegendreP[n,x],{n,1,10}],{x,0,1}]

Задача. Найдите производящую функцию для многочленов Лежандра.

T5=ChebyshevT[5,#]&;
T[x]
FullSimplify[T[Cos[x]]]

Задача. Найдите производящую функцию для многочленов Чебышева

T6 = ChebyshevT[6,#]&;
Integrate[T5[x] T6[x] (1-x^2)^(-1/2),{x,-1,1}]

Псевдослучайные величины

Histogram[RandomReal[NormalDistribution[0,1],200]]

ListPlot[{RandomReal[NormalDistribution[0,1],{500,2}],RandomReal[{-1,1},{200,2}]},Axes->None]

ListDensityPlot[RandomReal[NormalDistribution[0,1],{500,3}],Axes->None]

ArrayPlot[RandomReal[1,{30,40}]]

5. Программирование

Символические выражения

В системе Mathematica, все выражения имеют такую структуру: заголовок и несколько аргументов. В стандартной записи заголовок пишется вначале, а потом перечисляются аргументы через запятую внутри квадратных скобок
head[arg1,arg2,...]

FullForm[(x/(Sqrt[5]))+(y^2)+1/z]

TreeForm[%]

Head[%]

((x/(Sqrt[5]))+(y^2)+1/z)[[2]]

Clear[f]
x+y//f

3^(1/4)+1//N

{a,b,c}~Join~{d,e}

ReplacePart[a+b+c+d,3->x^2]

Подстановки

{x,x^2,a,b}/.x->3

f[a]+f[b]/.f[x_]->x^2

Списки

L={1,2,3}

L^2

L[[2]]

x^L
%/.x->5

Table[i^2,{i,6}]

Join[{a,b,c},{x,y},{u,v,w}]

Append[{1,2,3},4]

Sort[{xyz,3,2,a,1}]

Select[{1,2,4,7,6,2},EvenQ]

Select[{1,2,3,4,5,6},#>3&]

Select[{1,2,4,7,6,2},#>2&,1]

s[x_]=DivisorSigma[1,x]-x;
Select[Range[10000],s[#]==#&]

Задача. Среди всех натуральных чисел от 1 до 10000, выпишите все пары простых чисел - близнецов, то есть пары простых чисел вида p и p+2.

Задача. Среди всех натуральных чисел от 1 до 10000, выпишите все пары чисел x и y, таких, что s(x)=y и s(y)=x.

Length[{a,b,c,d}]

Задача Макар - Лиманова

op[l_]:=Select[{Length[l]}~Join~(l-1),#>0&]
Manipulate[Nest[op,{16},n],{n,1,30,1}]

Условия

abs[x_]:=If[x<0,-x,x]
abs/@{-1,0,1}

If[a<b,1,0]
If[a<b,1,0,Indeterminate]
If[TrueQ[a<b],1,0]

a=2;Which[a==1,x,a==2,b]

Циклы

For[i=0,i<4,i++,Print[i]]

For[i=1;t=x,i^2<10,i++,t=t^2+i;Print[t]]

For[i=1,i<1000,i++,If[i>10,Break[]]];i

Do[Print[i^2],{i,1,3}]

n=1;While[n<4,Print[n];n++]

t=x;Do[t=1/(1+k t),{k,2,6,2}];t

Nest[f,x,3]

Пример : заполненное множество Жюлиа многочлена z^2 - 1 - множество точек на комплексной плоскости, орбиты которых не убегают на бесконечность.

RegionPlot[Abs[Nest[#^2-1&,x+I y,20]]<2,{x,-1.7,1.7},{y,-.9,.9},AspectRatio->Automatic,MaxRecursion->5]

Модули

gcd[m0_,n0_]:=
Module[{m=m0,n=n0},
While[n≠0,{m,n}={n,Mod[m,n]}];
m
]

gcd[6,33]

fib[n_]:=
Module[{f},
f[1]=f[2]=1;
f[i_]:=f[i]=f[i-1]+f[i-2];
f[n]
]

fib[7]

Переменные, объявленные внутри модуля, являются локальными переменными.

t=7;
Module[{t},t=8;Print[t]]

t

Список упражнений, которые нужно сдать на занятиях

1. Проверьте формулу Лежандра для удвоения аргумента гамма - функции :

(Gamma[2 s])=(1/(Sqrt[π])) (2^(2 s-1)) (Gamma[s]) (Gamma[s+(1/2)])

2. Наберите эту формулу (чтобы она выглядела так, как показано выше).  
Указание: читайте Help/ Documentation Center/ Notebooks and Documents/ Math Typesetting.

3. Нарисуйте график функции

f(x, y) = sin(x + y) + cos(xy)

4. Для натуральных значений n, составьте таблицу явных значений ζ (-n) для n=1,...,20. Та же задача для ζ(2n).

5. Найдите (с точностью до 10 знаков после запятой) наименьшее число y>0, такое что ζ(1/2+yi)=0. 
Указание: при помощи команды FindRoot, найдите корень модуля дзета-функции (в качестве начальной точки укажите 0.1).

6.  У Вас в кармане 30 монет достоинством 1 рубль, 20 монет достоинством 2 рубля, 15 монет достоинством 5 рублей, и 10 бумажек достоинством 10 рублей. Булочка стоит 26 рублей. Сколькими способами Вы можете за нее расплатиться?
Указание: задача сводится к нахождению коэффициента при (x^26) в многочлене

p(x)=((∑_(k=0)^30)(x^k)) ((∑_(k=0)^20)(x^(2 k))) ((∑_(k=0)^15)(x^(5 k))) (∑_(k=0)^10)(x^(10 k))

7. Нарисуйте график логарифма числа Бернулли (Subscript[B, 2k]) в зависимости от натурального числа k.
Указание: читайте раздел Алгебра / Производящие функции выше.

8. Среди многочленов вида (x^n)+x+1, найдите неприводимые (над рациональными числами).
Указание: используйте функции Select и IrreduciblePolynomialQ.

9. Особенность "ласточкин хвост" возникает на поверхности, заданной в трехмерном пространстве с координатами a,b,c уравнением
	Discriminant[a (x^3)+b (x^2)+c x+(x^5),x]==0
Нарисуйте эту поверхность так, чтобы особенность в начале координат была хорошо видна.
Указание: используйте команду ContourPlot3D с подходящим значением параметра MaxRecursion

10. Гаусс обнаружил, что правильный 17-угольник можно построить циркулем и линейкой. Алгебраически, это означает, что sin((π/17)) представляется формулой, содержащей лишь квадратные радикалы. Выпишите эту формулу.
Указание: используйте команду ToRadicals.

11. Нарисуйте кривые уровня функции f∘(A^n), где f(x,y)=sin(x)+sin(y), a A=(({{1, 1}, {1, 0}})), при n=1,...,10.
Указание: используйте команды ContourPlot, MatrixPower.

12. Рассмотрим матрицу
A=(({{5, 4, 2, 1}, {0, 1, -1, -1}, {-1, -1, 3, 0}, {1, 1, -1, 2}})).
Найдите ее характеристический и минимальный многочлены. Приведите ее к жордановой нормальной форме.

13. Найдите инвариантное распределение вероятностей для случайного блуждания на графе
GraphPlot[{1->2,2->3,3->4,4->1,5->1}]

Указание. Воспользуйтесь функцией, определенной выше в Линейная алгебра / Собственные векторы и собственные значения / Цепи Маркова и случайные блуждания.

14. Найдите пятую производную функции sin(cos(x)).

15. Нарисуйте на одном графике линии уровня функции sin(x)+cos(y) и ее градиент. Сделайте так, чтобы при изменении функции в программе пришлось бы менять только ее определение. Указание: воспользуйтесь командами Show,  ContourPlot, VectorPlot.

16. Найдите длину эллипса с главными полуосями a, b. Указание: используйте Integrate с Assumptions -> {a > 0, b > 0, a < b}.

17. Найдите разложение котангенса в степенной ряд Лорана (до членов порядка 20).

18. Посчитайте кривизну и кручение кривой, заданной параметрически формулами x=t, y=(t^2), z=(t^3). Указание: см. 
http://en.wikipedia.org/wiki/Torsion_of_curves
Там есть формула, по которой можно посчитать кручение пространственной кривой.

19. Вычислите неопределенный интеграл∫(1/(sin(x)+cos(x)))dx

20.  Найдите (численно, с точностью до 4 знаков после запятой) площадь эллипсоида с главными полуосями 1, 2, 3. Указание: эллипсоид может быть запараметризован таким образом:
x = cos(u) cos(v)
y = 2 sin(u) cos(v)
z = 3 sin(v).

21. Найдите производящие функции для многочленов Лежандра и для многочленов Чебышева. Указание: есть команда GeneratingFunction.

22. Найдите число нулей многочлена (z^5)+5z-1 в круге (x^2)+(y^2)<1.2, z=x+iy.

23. Вычислите производную неявной функции y(x), заданной уравнением sin(xy)+(y^3)=0.
Производную нужно выразить через x и y.

24. Пусть f(n) - количество натуральных чисел x<n, таких, что числа (x^2)+1 простые. Определите функцию f.

25. Напишите функцию, представляющую четное натуральное число в виде суммы двух простых.

Проекты

Выпуклые оболочки

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

Число нулей дзета-функции

Посчитайте число нулей дзета-функции в прямоугольнике [0.25,0.75]x[0,M] для некоторых больших значений M (например, для M=10000). Оцените, как растет это число с ростом M.

Какие из многообразий Шуберта гладкие?

Всякой перестановке соответствует алгебраическое многообразие - многообразие Шуберта. Некоторые многообразия Шуберта гладкие, а некоторые - нет. По перестановке можно достаточно легко понять, гладко соответствующее многообразие Шуберта или нет. Алгоритм такой. Рассмотрим перестановку p. Если последовательность p(a)p(b)p(c)p(d) совпадает с последовательностью dbca или с  последовательностью cdab для некоторых a<b<c<d, то соответствующее многообразие Шуберта негладкое (а иначе гладкое).

Найдите все перестановки из четырех и пяти элементов, такие, что соответствующие многообразия Шуберта негладкие. Посчитайте число негладких многообразий Шуберта для перестановок из шести элементов.

Абелевы группы

Конечно порожденную абелеву группу можно задать образующими и соотношениями. Пусть, например, есть 10 образующих и 5 соотношений. Тогда все соотношения можно записать в виде целочисленной матрицы 10 на 5. При помощи элементарных целочисленных операций над строками и столбцами этой матрицы, ее можно привести к диагональному виду, причем так, что на диагонали будут стоять целые положительные числа, линейно упорядоченные отношением делимости. Эти числа (так называемые инвариантные множители) фигурируют в представлении нашей группы в виде прямой суммы циклических. Напишите программу, вычисляющую инвариантные множители для матриц 10 на 5.

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

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

Поверхности тора и кренделя

Нарисуйте поверхности тора и кренделя. Программа должна быть оформлена в виде манипулятора, позволяющего менять толщину и угол поворота.

Цепные дроби, ζ(3) и уравнение Пелля

Запрограммируйте манипулятор, который, в зависимости от целых чисел a,b и натурального числа n, выводил бы n первых элементов цепной дроби для числа

(π^a) (e^b) (Zeta[3]).

Цепная дробь должна выводиться в традиционной записи, а не в виде списка. Кроме того, напишите функцию, которая, по натуральному числу d, выдавала бы какое-нибудь решение уравнения Пелля (x^2)-d (y^2)=1 в целых числах (решение может быть получено при помощи цепных дробей, см. http://en.wikipedia.org/wiki/Pell_equation).

Бассейны притяжения для рациональных функций

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

(Subscript[f, a])(x)=(a/((x^2)+2 x))

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

Аппроксимации Лагранжа и Паде

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

Выпуклые многогранники

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

Вырождения критических точек

Критическая точка функции f называется вырожденной, если второй дифференциал в этой точке - вырожденная квадратичная форма. Для двух многочленов четвертой степени P и Q от двух переменных, найдите значения параметра t, при которых функция P+tQ имеет вырожденную критическую точку. Запрограммируйте манипулятор, который, по значению параметра t, рисует критические точки функции P+tQ, и указывает индексы этих критических точек (индекс = число минусов в диагонализации второй квадратичной формы). Та же задача для тригонометрических многочленов, определенных на торе.

Геометрическая динамика

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

Числа Гурвица

Напишите программу вычисления чисел Гурвица (Subscript[h, g,λ]), где g ≥ 0, а λ = ((Subscript[λ, 1])≥(Subscript[λ, 2])≥...≥(Subscript[λ, s]))  - разбиение. Это количество наборов из m=∑((Subscript[λ, i])+1)+2g-2 транспозиций в группе (Subscript[S, n]), n=(∑_( i))(Subscript[λ, i]), таких что их произведение имеет циклический тип λ. Количество
нужно разделить на n! Выпишите члены малых степеней в производящей функции для чисел Гурвица. Подробности см.лекции Ландо прошлого года по комбинаторике, есть на vyshka.math.ru

Ряды Эйзенштейна

Ряд Эйзенштейна определяется формулой

(Subscript[G, r])(τ)=(∑_(n=-∞)^∞)(∑_(m=-∞)^∞)(1/(((m τ+n))^r)),

в которой суммирование производится по всем парам целых чисел (n,m), таким что m и n не равны одновременно нулю. Число r - целое (ряд сходится при r>2). Очевидно, что для нечетных r ряды Эйзенштейна равны нулю. При четных r, ряд (Subscript[G, 2r])  выражается через ряды (Subscript[G, 4]) и (Subscript[G, 6]) согласно следующему рекуррентному соотношению:

(2 k+1) (k-3) (2 k-1) (Subscript[G, 2 k])=3 (∑_(r=2)^(k-2))(2 r-1) (2 k-2 r-1) (Subscript[G, 2 r]) (Subscript[G, 2 k-2 r])

Напишите программу, выражающую ряд (Subscript[G, 2r])  в виде многочлена от (Subscript[G, 4]) и (Subscript[G, 6]). При помощи этой программы, найдите (Subscript[G, 100]).

Замощения плоскости

Нарисуйте фундаментальную область группы замощений плоскости. Сделайте так, чтобы пользователь мог менять форму фундаментальной области (но чтобы она всегда оставалось фундаментальной областью), и чтобы мог выбирать одну из 17 групп замощений.

Синтаксический анализ

Напишите функцию, которая по текстовой строке, содержащей цифры, скобки и арифметические операции, производит указанное вичисление. Например, должно быть
eval["(2+3*5)-(6+10)"] = 1

Диаграммы Вороного

Рассмотрим конечное множество точек X на плоскости. Клетка Вороного элемента x∈X определяется как множество точек плоскости, расстояние от которых до x меньше, чем до любого другого элемента множества X. Диаграмма Вороного - это разбиение плоскости на клетки Вороного точек множества X. Напишите манипулятор, позволяющий ставить точки, или выбирать случайные наборы точек, и строящий по данным точкам диаграмму Вороного.

Векторные поля и предельные циклы

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

(dx/dt)=ϵ (-((x^3)/3)+x-y)
(dy/dt)=(x/ϵ)

Пользователь должен иметь возможность менять параметр ϵ.

Огибающие семейств прямых

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

Поризм Понселе

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

Идеалы в кольце многочленов

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