Педагогический вестник


Матвей ЭНТИН
Системы аналитических вычислений на примере программы "Mathematica"

При решении различных математических и физических задач для многих школьников и студентов большую проблему представляют громоздкие вычисления. Часто проблемой является не принципиальное понимание подхода к задаче, а именно трудоемкость вычислений. Учащийся запутывается в преобразованиях, делает арифметические ошибки. В результате задача оказывается нерешенной. Системы для аналитических вычислений представляют могущественные механизмы решения различных математических задач. Умение пользоваться ими не только помогает в обучении, но и дает мощный метод для профессиональной работы. Существует ряд таких систем, среди которых отметим MathCad, Derive, Reduce, Mapple и Mathematica. Автор имеет опыт использования последней системы в течение ряда лет и поэтому остановится именно на ней.

Системы, подобные программе "Mathematica", могут многое. Во-первых, они содержат арсенал обычных систем программирования типа C++ или Basic. Но к ним добавляется мощные графические системы, позволяющие строить графики в 2 и 3 измерениях. Имеется достаточно развитая численная математика, которая умеет решать численно алгебраические и дифференциальные уравнения. Но главным отличительным свойством является возможность проводить аналитические расчеты.

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

103678125216741057526102517377541798754344128431078736692246786/ 109179036977856376061845912489852645543544787726853619942498251/ 011107830946835203614112542972225659299697132701903790722407189/ 513569110221602223555117660629641350394018785824453797020371866/ 016015408069862400386359555153564259844545538088763279906102068/ 292419906072015115957670323764805727044793901151113811651794598/ 819565862871157884595738561896854046434800778773574599234968869/ 087824889020339785550030018212709889337497147584853266659939453/ 669500880268248919650797492728660855460141759701301568069771121/ 049327141010079690488477116488026912384231430471176350504501116/ 525908227572596249241035975472963490584534466543182183393095093/ 233233575270747221439688510833338211771660996008791411908629500/ 108351674101463810522158684028092503318669113053769277323826813/ 528766221557559559173068693279944870867705265245882873058130307/ 108395029065056431837398674449

\ - знак переноса на другую строку.

Теперь арифметика. Упростим дробь17017/20111. Если попытаться сделать это численно, то получим 0.846154. А Mathematica способна упростить это выражение: Simplify[17017/20111] = 11/13. Мы применили функцию "упростить"- Simplify[] к дроби. Стандартные функции в Mathematica начинаются с большой буквы и аргумент заключается в квадратные скобки.

Теперь упростим алгебраическое выражение
a4+4a3b+6a2b2+4ab3+b4:

Simplify[a^4 + 4a^3 b + 6 a^2 b^2 + 4 a b^3 + b^4] = (a+b)^4.

Значок ^ обозначает степень: a^b - это a в степени b.

Теперь решим алгебраическое уравнение. Это делает функция Solve (решить). Функция Solve имеет вид Solve[a==b,x]. Первый аргумент функции - это само уравнение. Двойной знак равенства - не ошибка, а обозначение уравнения. Второй аргумент - переменная, относительно которой уравнение a==b следует решать.

Решаем уравнение x1/2 + x - 1 = 0:

Solve[x^(1/2)+x-1==0,x]

Результат: (3-5^(1/2))/2.

Другой пример - тригонометрическое уравнение sinx + cosx = 1:

Solve[Sin[x]+Cos[x]==1,x]

Результат: {{x ->0}, {x ->Pi/2}}. Фигурные скобки обозначают массив значений, а стрелки, какие значения принимают решения. Pi - это число p. Математика, правда, предупредила, что найденные решения - не все из имеющихся. Для проверки построим график функции sinx + cosx - 1:

Plot[Sin[x]+Cos[x]-1,{x,0,5Pi}]

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



Видна периодичность функции 2p, свидетельствующая о том, что решение следует искать в виде 2pn или p/2+2pn (это только наводящее соображение, требуется корректное доказательство).

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

x + y + z = 1,
x - 2y + z = 0,
2x + y = 3.

Solve[{x+y+z==1,x-2y+z==0,2x+y==3},{x,y,z}]

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

Результат: {{x -> 4/3, y->1/3, z ->-2/3}}.

Физический пример. Пусть мы хотим найти сопротивление цепочки сопротивлений R1 = 1, R2 = 2, R3 = 1, R4 = 1, R5 = 3 Ом.



Составляем уравнения из законов Кирхгофа:

U = U1 + U2 = U4 + U5, I = I1 + I4, I1 = I2 + I3, I4 + I3 = I5, U1 + U3 = U2, Uj = RjIj.

Или

U = I1 + 2I2 = I4 + 3I5, I = I1 + I4, I1 = I2 + I3, I4 + I3 = I5, I1 + I3 = I4.

Решаем уравнения:

Solve[{i1+2i2==R,i4+3i5==R,i1+i4==1,i1==i2+i3,i4+i3==i5,i1+i3==i4},{i1,i2,i3,i4,i5,R}]

Результат {{i1 -> 9/17, i2 -> 10/17, i3 -> -1/17, i4 -> 8/17, i5 -> 7/17, R -> 29/17}}

Рассмотрим численное решение дифференциальных уравнений механики с помощью системы Mathematica. В качестве примера решим уравнение Ньютона для колебаний осциллятора:

p'(t) = -kx(t),
mx'(t) = p(t).

Здесь k - жесткость пружины, m - масса частицы, x(t) - отклонение в момент времени t, p(t)-импульс. Для решения уравнения его нужно дополнить начальным условием. Пусть в начальный момент времени частица отклонена до точки 1: x(0) = 1, а скорость ее равна нулю: p(0) = 0. Примем массу и жесткость за единицу. Вся программа на Mathematica для построения траектории колебаний имеет вид:

m=1;k=1;
eq={p'[t]==-k x[t], m x'[t]==p[t],x[0]==1,p[0]==0}
nds=NDSolve[eq,{x[t],p[t]},{t,0,100}];
Plot[x[t]/.nds,{t,0,100},PlotRange->All]

Первая строчка конкретизирует значения m и k. Во второй строчке сформулирована система уравнений и начальных условий, обозначенная eq. Эта система - список, состоящий из отдельных уравнений. Элементы списка в Mathematic'е разделеляются запятыми и заключаются в фигурные скобки. Третья строчка - функция NDSolve, находящая решение дифференциального уравнения eq относительно функций {x[t],p[t]} на промежутке времени t от 0 до 100. Три аргумента функции NDSolve разделены запятыми. Четвертая строчка - построение графика функции x[t], найденной при помощи решения уравнения nds (это выражено с помощью операции /.nds). График нужно построить на промежутке 0 < t < 100. Это выражено списком во втором аргументе функции Plot. Третий (необязательный) аргумент указывает, что на график должны попасть все точки кривой. Результат - график колебаний.



Модифицируем первое уравнение, введя жидкое трение (сила трения пропорциональна скорости Fтр = -bv):

eq={p'[t]==-k x[t]-(b /m)p[t], m x'[t]==p[t],x[0]==1,p[0]==0}.

Колебания становятся затухающими.




Теперь нарисуем зависимость импульса от координаты. Для этого воспользуемся параметрическим графиком ParametricPlot:

ParametricPlot[{x[t],p[t]}/.ds,{t,0,100},PlotRange->All,AspectRatio->1]

Первый аргумент - список двух координат, второй аргумент - список-параметр t и его интервал изменения 0,100. Четвертый необязательный аргумент определяет соотношение масштабов по осям y и x. Без трения график - окружность.




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




Рассмотрим еще один пример: стреляем из пушки. Полет снаряда описывается уравнением Ньютона. Обычно в школе разбирается движение брошенного камня в отсутствии сопротивления воздуха. Траектория снаряда в этом случае является параболой и описывается простой формулой из учебника. А как выглядит траектория полета с учетом сопротивления воздуха? К сожалению, решить уравнение движения снаряда в этом случае аналитически сложно. А промоделировать на компьютере очень просто.

Мы заложим модель сопротивления воздуха, исходя из того, что на снаряд действует сила давления набегающего потока воздуха arv2S = bv2, где r - плотность воздуха, v - скорость снаряда, S - поперечное сечение снаряда, a - численный коэффициент, связанный с формой снаряда. Сила сопротивления направлена против движения снаряда. Вся программа на Математике имеет вид

m=1;g=10;b=0.01;p=200;f=1;p1=p Cos[f];p2=p Sin[f];
eq={m x'[t]==px[t],m y'[t]==py[t],px'[t]==-b px[t](px[t]^2+py[t]^2)^(1/2)/m,py'[t]==-m g-b py[t](px[t]^2+py[t]^2)^(1/2)/m, x[0]==0,y[0]==0,px[0]==p1,py[0]==p2};
ns=NDSolve[eq,{x[t],y[t],px[t],py[t]},{t,0,100}];
ParametricPlot[{x[t],y[t]}/.ns,{t,0,20},PlotRange->All,AspectRatio->Automatic]

В первой строчке определены параметры: масса снаряда m, ускорение силы тяжести g, начальный импульс p, константа g, начальный угол к горизонтали f (в радианах), начальные горизонтальный и вертикальный импульсы p1 и p2. Система уравнений Ньютона с начальными условиями дается выражением eq; px[t] и py[t] - горизонтальный и вертикальный импульсы; ns - это функция, решающая систему уравнений Ньютона eq. Ниже приведены типичные траектории полета. В отличие от случая с отсутствием сопротивления воздуха, траектория несимметрична: снаряд, затормозив поступательную скорость, начинает круто падать вниз.




Некоторые тонкости Mathematic'и.

Синтаксис в Mathematic'е свой. Самое необычное - это то, что она помнит все присвоения, сделанные во время сессии. Поэтому не забывайте "освобождать" переменные перед вычислениями. Это делается присвоением типа x=. (x равно точке). Операторы можно отделять друг от друга точкой с запятой (;). Этот знак, кроме прочего, подавляет вывод на экран результата действия (кроме графики!) Функции определяются с неопределенным аргументом, обозначаемым аргументом с подчерком: F[x_]= x^2. Это означает, что в качестве аргумента может выступать любая величина.

Арифметика. Как обычно, работают знаки + ,-, * (умножить), /. Возведение в степень обозначается знаком ^: a^b значит a в степени b. Вместо знака * достаточно сделать пропуск или несколько. Но две буквы рядом будут восприняты как единый идентификатор. Цифра перед буквой воспринимается как умножение. Перед скобкой также необязателен знак умножения. Знак равенства = это присвоение. Выражение в правой части равенства вычисляется и присваивается выражению в правой части. Знак := обозначает задержанное присвоение, в котором правая часть не вычисляется.

Численные и приближенные вычисления осуществляются с помощью функции N[выражение]. Функция N[выражение, n] дает приближенное значение с точностью до n разрядов. Если этой функции нет, Mathematic'а вычисляет выражение точно.

Выделить числитель и знаменатель дроби можно функциями Numerator[a/b]=a или Denominator[a/b]=b. Эти функции работают как над числовыми дробями, так и над алгебраическими дробями.

Функция Together[aa] приводит выражение к общему знаменателю. Так, Together[a/b+c/d]=(ad+bc)/bd. Функция Expand[aa] разлагает выражение на простые дроби: Expand[(a d+b c)/bd]= a/b+c/d.

Помощь в Mathematic'е. В современных версиях имеется хорошо развитая система help'ов. Имеется книга создателя системы Вольфрама. Имеются демонстрации, поиск по ключевым словам. Из help'а можно копировать примеры в рабочее окно, а также изменять программу и вычислять прямо в help-окне. Mathematica Tour (тур) позволяет получить общее представление о возможностях системы.

Я привел только крохотную часть возможностей Mathematica. Более старые версии Mathematica работают под Windows 3.0 и даже под DOS, для работы с которыми могут подойти и 386 и 286 компьютеры.

Я рекомендую всем изучающим физику освоить эту систему. Это очень полезно для работы!


[ Предыдущее сообщение      Оглавление     Последующее сообщение]



vlad@ssl.nsu.ru