Мостовской А.П. Информационные технологии в математике

Подождите немного. Документ загружается.

1.24 Программирование в Математике 71

EllipticTheta[i, z, q] - эллиптическая θ-функция: θ

(z, q), i = 1, 2,

3, 4.

EulerE[n, x] - n-й полином Эйлера.

EulerPhi[n] - Эйлерова функция ϕ(n) ( количество положительных

целых чисел, не превосходящих n и взаимно простых с n).

ExpIntegralE[n, z] - интегральная показательная функция E(n,z).

ExpIntegralEi[z] - интегральная показательная функция Ei(z).

LegendreP[n, m, x] - присоединенный полином Лежандра.

LegendreQ[n, z] - n-ая функция Лежандра I рода.

LegendreQ[n, m, z] - присоединенная функция Лежандра II рода.

CoshIntegral[x] - гиперболический интегральный косинус.

LogIntegral[z] - интегральный логарифм li(z).

SinhIntegral[x] - интегральный гиперболический синус

Integrate[Sinh[t]/t, {t, 0, x}].

SinIntegral[x] - интегральный синус Integrate[Sin[t]/t, {t, 0, x}].

Gamma[a] - Эйлерова гамма-функция Γ(a).

Hypergeometric0F1[a, z] - гипергеометрическая функция 0F1(a, z).

JacobiSymbol[n, m] - символ Якоби от n и m.

JacobiZeta[ϕ, m] - дзета-функция Якоби Z(ϕ|m).

MoebiusMu[n] - функция Мебиуса µ(n).

SphericalHarmonicY[l, m, θ, ϕ] - сферическая гармоника Ylm(θ, ϕ).

StirlingS1[n, m] - число Стирлинга I рода.

StirlingS2[n, m] - число Стирлинга II рода.

WeierstrassP[u, g2, g3] - эллиптическая функция Вейерштрасса P.

Zeta[s] - дзета-функция Римана ζ(s).

1.24 Программирование в Математике

Математику можно рассматривать как язык программирования вы -

сокого уровня, ориентированный на решение математических задач.

Все команды, рассмотренные выше, являются и командами этого языка

программирования. Рассмотрим базовые команды программирования

- команды цикла и команды ветвления в Математике.

1.24.1 Циклическое выполнение команд

Циклы в Математике организуется с помощью команды Do:

72 ГЛАВА 1. СИСТЕМА МАТЕМАТИКИ

Do[expr, {m}] выполняет m раз команду expr. Здесь и далее под

expr понимается команда или несколько команд, разделенных точкой

с запятой.

Пример: s = 0; Do[s = s + 1, {100}]; s ⇒ 100. Здесь три команды

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

Do[expr, {i, m}] выполняет expr m раз, при этом целая переменная

i принимает значения от 1 до m (с шагом 1).

Например, найдем сумму первых 100 натуральных чисел:

s = 0; Do[s = i + s, {i, 100}]; s ⇒ 5050.

Do[expr, {i, a, b}] - выполняет expr всякий раз, как только изме-

няется счетчик циклов i от a до b с шагом 1.

Пример суммирования по i от 12.2 до 23.5 с шагом 1:

s = 0; Do[s = i + s; Print[i], {i, 12.2, 23.5}]; s ⇒ 212.4

Do[expr, {i, a, b, d}] - счетчик циклов i принимает значения от a

до b с шагом d, всякий раз выполняется команда expr.

Do[expr, {i, i1, i2}, {j, j1, j2}, ...] - организуется цикл по j для

каждого i и так далее.

Пример. Следующий цикл выполняется по двум индексам i,j:

s = 0; Do[s = j+i+s; P rint[”i = ”, i, ””, ”j = ”, j], {i, 1, 3}, {j, 1, 4}]; s

Другой вид цикла реализуется командой For:

For[start, test, incr, body] выполняет команду start, затем ци-

клически вычисляет body и incr до тех пор, пока test не перестанет

давать логическое значение True.

Например, F or[s = 1, s < 6, s + +, P rint[s]] или

F or[i = 1; t = x, i

< 10, i + +, t = t

+ i; P rint[Expand[t]]].

Здесь s++ ( или s+=1) означает, что s=s+1. Аналогично, s– ( или

s-=1) означает, что s=s-1.

Цикл можно также организовать с помощью функции While:

While[test, expr] выполняет expr до тех пор, пока test дает логи-

ческое значение True.

x=1; W hile[x < 30, x = x + 1; y = x

; P rint[y]].

1.24.2 Команды ветвления

Для организации разветвляющихся программ применяются следую-

щие команды:

If[condition, t, f]

возвращает t, если результатом вычисления condition будет True,

1.24 Программирование в Математике 73

и возвращает f, если результат False.

If[condition, t, f, u]

возвращает t, если результатом вычисления condition будет True,

возвращает f, если результат False и возвращает u, если в результате

вычисления condition не будет получено ни True, ни False.

Примеры.

1. If [1 > 2, a, b] ⇒ b.

2. If [x == y, a, b] ⇒ If[x == y, a, b]. Если x,y не определены, то

уравнение x==y не есть высказывание, поэтому Математика непонят-

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

3.If[x == y, a, b, c] ⇒ c.

4. Определим функцию f[x ] := If[x > 0, x

, −x

] и построим ее

график Plot[f[x],{x,-2,2}].

5. Напомним, что тройное равенство p===q означает, что p и q

синтаксически тождественны. Определим функцию

g[x ] := If[Head[x] === P lus, x

, x

]

и найдем значения g[a+b]; g[a]. Аргументом функции g может быть

любое выражение. Если заголовок Head[x] выражения x есть Plus (это

так, если x=a+b), то g[a + b] = (a + b)

, а если заголовок x не Plus,

(Head[a]=Symbol), то g[a] = a

Множественное ветвление может быть организовано с помощью ко-

манд Which и Switch:

Which[test

, value

, test

, value

, ...] - команда проверяет, начиная

с первого, все высказывания test

. Если i-ое высказывание test

впер-

вые приняло значение True, то команда возвращает value

Switch[expr, form

, value

, form

, value

, ...] вычисляет expr, срав-

нивает его последовательно с каждым form

, затем вычисляет и воз-

вращает то value

, которое соответствует первому совпадению.

Примеры.

1. W hich[1 == 2, 1, 2 == 2, 2, 3 == 3, 3] ⇒ 2.

2. Определим функцию - модуль числа f[x] = |x| : f[x ] := W hich[x <

0, −x, x >= 0, x]. Построим ее график: Plot[f[x],{x,-3,3}]. Эти две коман-

ды можно объединить в одну:

P lot[W hich[x < 0, −x, x >= 0, x], {x, −3, 3}].

3. Определим функцию: f[x

] := Swtch[x, Integer, x

, Real, x

]. Фун-

кция проверяет тип числа x: f[3]=27,f[5.5]=30.25.

Задачи. 1). Функцию f[x ], определенную на промежутке [0,4] ра-

венствами:

74 ГЛАВА 1. СИСТЕМА МАТЕМАТИКИ

f[x ] := x(2 − x)/; ((x >= 0)&&(x <= 1))

f[x ] := 1/; ((x >= 1)&&(x <= 2))

f[x ] := (4 − x)/2/; ((x >= 2)&&(x <= 4))

задать с помощью команды Which и найти ее производную, первообраз-

ную и определенный интеграл.

2). Найти n-ое число Фибоначчи.

1.25 Анимация и звук

Математика позволяет проследить за изменением геометрического

образа в реальном режиме времени. Пусть, например, график pic[t]

зависит от параметра t. Команда

Do[pic[t],{t,a,b,d}]

воспроизводит на экране последовательно pic[a], pic[a+d],...,. Если

теперь дважды щелкнуть мышкой по одному и тому же графику или

мышкой затемнить скобку в правой части экрана, содержащую все эти

графики (то есть щелкнуть правой клавишей мышки , когда ее указа-

тель установлен на скобке, содержащей все графики) и нажать последо-

вательно клавиши Ctr+Y, то графики на экране начнут самостоятель-

но последовательно меняться. Частота кадров управляется мышкой и

клавишами небольшой панели, появившейся внизу экрана. Чтобы от-

менить анимацию, достаточно щелкнуть мышкой в рабочем поле . В

Математике 4 для запуска анимации надо щелкнуть дважды мышкой

на одном из графиков.

Примеры.

1. Do[ContourP lot[Sin[x

+ i] − Cos[y

+ i], {x, −3, 3}, {y, −3, 3},

ColorF unction−> Hue, P lotP oints−> 525, F rame−> F alse], {i, 1, 10}]

2. Do[P lot3D[Cos[Sqrt[x

+ y

] + Abs[t − 2π]]/Sqrt[x

+ y

+ 1/4]

{x, −4π, 4π}, {y, −4π, 4π}, P lotP oints−> 26, P lotRange−> {−2, 2},

BoxRatios−> {1, 1, 1}, Boxed−> F alse, Axes−> F alse],

{t, 0, 2π − (2π/16), 2π/16}].

Математика имеет возможность синтезировать звук. Звук можно

воспроизвести, если компьютер оборудован звуковой картой. Возможен

синтез как монофонических, так и стереофонических звуков. С синте-

зируемым звуком связан некоторый графический образ - ячейка. Он

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

Для синтеза звуков используются следующие команды:

ListPlay[{a1, a2, ...}]

1.26 Две задачи математической физики 75

проигрывает звук с амплитудой, заданной последовательностью уров-

ней ai,

Play[f[t], {t, a, b}]

воспроизводит звук с амплитудой, заданной как функцией f[t] от вре-

мени t ∈ [a, b].

Команда

Play[{f[t],g[t]} {t, a, b}]

воспроизводит стерео звук. Левый канал описывается функцией f[t].

Пример. P lay[Sin[

1000

], {t, −2, 2}].

1.26 Две задачи математической физики

1.26.1 Колебание длинной струны. Метод Даламбера

Пусть функция

f[x

] := Which[x <= −2, 0, −2 < x < 0, (1/4)x + 1/2,

0 < x < 2, −x/4 + 1/2, x > 2, 0]

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

времени t=0. (Можно посмотреть на график этой функции:

P lot[f[x], {x, −20, 20]}).

Скорость g[x] точек струны при t=0 будем считать равной нулю

- g[x ] := 0. Известно, что колебание струны под действием сил на-

тяжения и сил, связанных с ускорением g[x], описывается формулой

Даламбера:

u[x , t ] :=

(f[x − at] + f[x + at])/2 + (1/(2a))Integrate[g[z], {z, x − at, x + at}],

где a - параметр струны.

Рассмотрим колебание струны как анимацию решения с заданной

скоростью:

a = 3; Do[Plot[u[x, t], {x, −20, 20}, Ticks−> {Automatic, {0, 0.2}},

AspectRatio−> 1/3]; Print[”t = ”, t], {t, 0, 3, 0.05}].

1.26.2 Колебание струны с закрепленными концами

Пусть функция f(x), определенная на интервале [0,l], за дает положе-

ние плоской струны в начальный момент времени t=0. Пусть концы

струны закреплены в точках x = 0 и x = l оси x-ов.

Колебание такой струны описывается функцией u[x,t], равной сумме

76 ГЛАВА 1. СИСТЕМА МАТЕМАТИКИ

бесконечного ряда Фурье:

u(x, t) =

∞

k=1

cos

t + b

sin

t) sin

При этом граничные условия задаются равенствами u(0, t) = u(l, t) =

0, начальные условия - u(x, 0) = f(x), u

(x, 0) = g(x), а коэффициенты

ряда вычисляются по формулам:

f(x) sin (

x)dx, b

akπ

g(x) sin (

x)dx

Рассмотри конкретную задачу. Введем команду Clear[f,g,u,a].

Пусть функция f[x

] := Which[x <= 4, (1/4)x, x > 4, −x + 5], опре-

деленная на интервале [0,5], задает положение плоской струны в на-

чальный момент времени t=0. (Можно посмотреть на график этой

функции:

P lot[f[x], {x, 0, 5]}).

Пусть скорость точек струны при t=0 равно нулю: g[x ] := 0.

Вычислим коэффициенты ряда:

a[k ] := (2/5)Integrate[(1/4)x Sin[k π x/5], {x, 0, 4}]+

(2/5)Integrate[(−x + 5)Sin[k π x/5], {x, 4, 5}];

b[k ] := 0.

Возьмем функцию u[x,t] приближенно - в виде суммы первых 4 чле-

нов ряда:

u[x , t ] = Sum[a[k]Cos[t a k π/5]Sin[k π x/5], {k, 1, 4}].

Зададим анимацию :

a = 1; Do[Plot[u[x, t], {x, 0, 5}, Ticks−> {Automatic, {0, 0.2}},

AspectRatio−> 1/3], {t, 0, 20, 0.4}].

Перед повторным вводом приведенных команд или при изменении

параметра струны a надо ввести коман ду очищения: Clear[f,g,u,a].

Задачи. 1). Рассмотреть анимацию решения дифференциального

уравнения, описывающего колебание бесконечной струны, с начальны-

ми условиями: a) u(х,0)=x, u

(x, 0) = −x, b) u(х,0)=0, u

(x, 0) = cos(x).

2). Рассмотреть анимацию решения дифференциального уравнения

колебания струны, если a) струна закреплена на концах x=0, x=3. В на-

чальный момент времени струна имеет вид ломанной OAB, где O(0,0),

A(2,-0.1), B(3,0). Начальная скорость равна нулю.

b) струна закреплена на концах x=0, x=2. В начальный момент вре-

мени струна имеет вид отрезка [0,2], а начальное ускорение задается

функцией u

(x, 0) = cos

π(x−1)

Глава 2

Программирование в Математике

Программирование математических задач одна из замечательных

возможностей системы Математика. Все коман ды Математики вместе

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

рования - ветвления, цикла, образуют удобную программную среду, в

которой можно решить большинство математических задач символь-

ного и численного программирования .

2.1 Примеры программирования в Математике

Приведем несколько примеров программирования в Математике.

2.1.1 Перестановка строк матрицы

Поменяем местами в данной матрице m строки с номерами i и j. По-

строим функцию с аргументами - матрица m и номера строк, которые

следует переставить:

r[m

, i , j ] := Block[{p = m}, p[[i]] = m[[j]]; p[[j]] = m[[i]]; p]

Для организации функции применена команда Block, выделяющая

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

ная p, ей присвоено значение аргумента команды. Напомним, что в

теле Block нельзя ничего делать с аргументом, он там является защи-

щенной переменной, но с копией аргумента (p) можно делать все.

Пример.

Пусть m = {{1, 2, 3}, {4, 5, 6}}. Тогда r[m, 1, 2] ⇒ {{4, 5, 6}, {1, 2, 3}}.

2.1.2 Транспонирование матрицы

Приведем программу транспонирования квадратной матрицы m:

78 ГЛАВА 2. ПРОГРАММИРОВАНИЕ В МАТЕМАТИКЕ

t[m ] := Block[{p = m}, Do[p[[i, j]] = m[[j, i]], {i, 1, Length[m]},

{j, 1, Last[Dimension[m]]}]; p}]

Отметим, что команда Transpose[m], встроенная в ядро, также транс-

понирует матрицу.

2.1.3 Перестановка столбцов матрицы

Команды перестановки строк матрицы и транспонирования матри-

цы позволяют менять местами столбцы матрицы:

c[m

, i , j ] := t[r[t[m], i, j]]

Эту команду можно применять и в следующем виде:

c[m , i , j ] := Transpose[r[Transpose[m], i, j]]

2.1.4 Решение системы линейных уравнений методом Гаусса

Следующая команда, дополняющая три предыдущие (их надо рас-

сматривать, как одну программу), позволяет решить систему линей-

ных алгебраических уравнений методом Гаусса с выделением главного

элемента и выводом на экран промежуточных шагов вычислений. Про-

грамма находит решение системы (любого числа переменных), име-

ющей единственное решение и квадратную матрицу . При работе в

Mathematica 2.0 надо убрать из программы ”//TraditionalForm”. Про-

грамма состоит из выделенных нумерованных строк, комментарии не

принадлежат программе (полный текст программы без номеров и ком-

ментариев - в следующем параграфе).

Пусть w - матрица из коэффициентов уравнений и свободных чле-

нов, размерность матрицы ld × (ld + 1), где ld - число неизвестных и

уравнений системы. Определим основную функцию common аргумента

w. Она состоит из команд строк с 1 по 28, ограниченных фигурными

скобками и разделенных точками с запятыми.

1. common[w

] := {Clear[a1]; w1 = w;

в строках 2-6 из матрицы w выделяются матрица системы m и столбец

свободных элементов b, определяются ранг ngw расширенной матрицы

w и ранг ngmo матрицы m. Команда Dimensions[w] возвращает размер-

ность матрицы w в виде списка {ld, ld + 1}, команда Last выделяет из

списка последний элемент ld+1, из которого вычитается размерность

нулевого пространства матрицы. В результате получаем ранг матри-

цы w:

2. ngw = Last[Dimensions[w]] − Length[NullSpace[w]];

2.1 Примеры программирования в Математике 79

транспонируем матрицу w и выделяем у нее последнюю строку, полу-

чаем список, состоящий из свободных членов системы:

3. b = Last[Transpose[w]];

отбросив строку из свободных членов системы, получим матрицу си-

стемы m:

4. m = Transpose[Drop[Transpose[w], −1]];

5. ld = Length[w1];

6. ngmo = Last[Dimensions[m]] − Length[NullSpace[m]];

команды следующих пяти строк проверяют систему на совместность

(команда Abort[] прекращает выполнение программы):

7. If[ngw == ngmo == ld, Print[””],

8. If[ngw == ngmo < ld,

Print[”Бесконечно много решений.”]; Abort[],

Print[”Нет решений”]];

Abort[]];

введем список a1, элементы которого есть индексы при неизвестных

системы:

9. a1 = Array[a, {ld}];

строка 10 - начало основного цикла, конец цикла - строка 21, счет-

чик цикла d. Рассмотрим первый проход цикла - считаем, что везде в

строках 10-20 d=1.

10. Do[l = ld − d + 1;

в строке 1 введена копия w1 матрицы w, с которой и будем иметь де-

ло (Математика запрещает изменять матрицу-аргумент w). Находим

максимальный модуль коэффициентов при неизвестных:

11. q = Table[Abs[w1[[i + d − 1, j + d − 1]]], {i, 1, l}, {j, 1, l}];

12. s = Max[q];

и выясняем позицию (j5,i5) максимального элемента s в матрице w1, у

которой отбрасываем столбец свободных коэффициентов:

13. i5 = Position[Abs[Transpose[Drop[

Transpose[w1], ld + 1]]], s][[1, 2]];

14. j5 = Position[Abs[Transpose[Drop[

Transpose[w1], ld + 1]]], s][[1, 1]];

переставляем столбцы и строки матрицы w1 так, чтобы максимальный

элемент занял верхнее левое положение :

15. w1 = r[w1, d, j5];

16. w1 = c[w1, d, i5];

в списке номеров неизвестных системы a1 произведем перестановку,

80 ГЛАВА 2. ПРОГРАММИРОВАНИЕ В МАТЕМАТИКЕ

соответствующую перестановке столбцов:

17. a1 = r[a1, d, i5];

делим элементы первой строки матрицы w1 на максимальный элемент

и заменяем первую строку в матрице w1 на измененную:

18. w1 = w1/.w1[[d]]−> w1[[d]]/w1[[d, d]];

следующий цикл (строка 20) обнуляет все элементы первого столбца

матрицы w1, кроме первого, который равен 1. Полученная измененная

расширенная матрица системы w1 выводится на экран:

19. Do[

w1 = w1/.w1[[kk]]−> (w1[[kk]] − w1[[d]]w1[[kk, d]]),

{kk, d + 1, ld}];

20. Print[d, ” Шаг : ”, w1//TraditionalForm],

первый проход цикла закончился. Далее счетчику d присваивается зна-

чение 2 (строка 21) и все команды в строках 10-20 выполняются при

этом значении d. При этом в измененной матрице w1 будут рассматри-

ваться все строки и столбцы d по ld-1 включительно. Цикл проработает

ld-1 раз.

21. {d, 1, ld − 1}];

В следующей строке делим строку ld матрицы w1 на последний коэф-

фициент w1[[ld,ld]]:

22. w1 = w1/.w1[[ld]]−> w1[[ld]]/w1[[ld, ld]];

выводим на экран матрицу w1, приведенную к треугольному виду:

23. Print[”Ответ : ”, w1//TraditionalForm];

и начинаем находить решения системы. Через y[i] обозначены решения

системы с треугольной матрицей, то есть y[1] - первое неизвестное,...,

y[ld] - последнее, его сразу можно найти:

24. y[ld] = w1[[ld, ld + 1]]/w1[[ld, ld]];

находим остальные решения:

25. Do[y[i] = w1[[i, ld + 1]]−

Sum[w1[[i, j]]y[j], {j, ld, i + 1, −1}], {i, ld − 1, 1, −1}];

пусть список xo состоит из решений x[i] первоначальной системы. Ре-

шения x и y отличаются только порядком расположения элементов:

26. xo = Array[x, ld];

число p=a1[[1]][[1]] равняется номеру переменной x[1], стоящей на пер-

вом месте в a1, то есть x[p]=y[1] и так далее. Реализуем эти присваи-

вания следующим циклом:

27. Do[x[a1[[i]][[1]]] = y[i], {i, 1, ld, 1}];

и напечатаем окончательный результат: