Решение физических и математических задач с Maxima

Доступная литература и сеть Интернет в качестве "электронного помощника" студентов и школьников обычно позиционирует пакет MathCad, изредка — Maple или Mathematica. Материал данной главы содержит ряд разнородных задач, которые решались разными авторами вручную или при помощи MathCad.

7.1 Операции с полиномами и рациональными функциями

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

7.1.1 Упрощение алгебраических выражений

Пример:

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

(%i1)	g:(1/a-1/(b+c))/(1/a+1/(b+c))*(1+(b^2+c^2-a^2)/2/b/c)/((a-b-c)/a/b/c);

$$(\%o1)\quad \frac{a\,b\,c\,\left( \frac{1}{a}-\frac{1}{c+b}\right) \,\left( \frac{{c}^{2}+{b}^{2}-{a}^{2}}{2\,b\,c}+1\right) }{\left( -c-b+a\right) \,\left( \frac{1}{c+b}+\frac{1}{a}\right) }$$$

(%i2)	ratsimp(%);

$$(\%o2)\quad -\frac{a\,c+a\,b-{a}^{2}}{2}$$

(%i3)	%,a=0.02,b=-11.05,c=1.07;

$$(\%o3)\quad 0.1$$

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

(%i1)	(sqrt(x)+1)/(x*sqrt(x)+x+sqrt(x))/(1/(x^2-sqrt(x)));

$$(\%o1)\quad \frac{\left( \sqrt{x}+1\right) \,\left( {x}^{2}-\sqrt{x}\right) }{{x}^{\frac{3}{2}}+x+\sqrt{x}}$$

(%i2)	ratsimp(%);

$$(\%o2)\quad x-1$$

Пример: Сделать указанную подстановку и результат упростить:

(%i3)	expr:(x^3-a^(-2/3)*b^(-1)*(a^2+b^2)*x+b^(1/2))/(b^(3/2)*x^2);

$$(\%o3)\quad \frac{{x}^{3}-\frac{\left( {b}^{2}+{a}^{2}\right) \,x}{{a}^{\frac{2}{3}}\,b}+\sqrt{b}}{{b}^{\frac{3}{2}}\,{x}^{2}}$$

(%i4)	ratsimp(%);

$$(\%o4)\quad \frac{{a}^{\frac{2}{3}}\,b\,{x}^{3}+\left( -{b}^{2}-{a}^{2}\right) \,x+{a}^{\frac{2}{3}}\,{b}^{\frac{3}{2}}}{{a}^{\frac{2}{3}}\,{b}^{\frac{5}{2}}\,{x}^{2}}$$

(%i5)	radcan(%);

$$(\%o5)\quad \frac{{a}^{\frac{2}{3}}\,b\,{x}^{3}+\left( -{b}^{2}-{a}^{2}\right) \,x+{a}^{\frac{2}{3}}\,{b}^{\frac{3}{2}}}{{a}^{\frac{2}{3}}\,{b}^{\frac{5}{2}}\,{x}^{2}}$$

Без указанной подстановки упрощение посредством комбинации функций $$ratsimp$$ и $$radcan$$ не удаётся.

(%i6)	%,x=a^(2/3)*b^(-1/2);

$$(\%o6)\quad \frac{\frac{{a}^{\frac{2}{3}}\,\left( -{b}^{2}-{a}^{2}\right) }{\sqrt{b}}+{a}^{\frac{2}{3}}\,{b}^{\frac{3}{2}}+\frac{{a}^{\frac{8}{3}}}{\sqrt{b}}}{{a}^{2}\,{b}^{\frac{3}{2}}}$$

(%i7)	ratsimp(%);

Конечный результат оказывается простым

$$(\%o7)\quad 0$$

7.1.2 Разложение полиномов и рациональных выражений на множители

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

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

Пример (также из сборника под ред. М.И. Сканави): Решить уравнение $$\sqrt(x-2) = x-4$$:

(%i1)	solve([sqrt(x-2)=x-4],[x]);

$$(\%o1)\quad [x=\sqrt{x-2}+4]$$

Уравнение имеет одно решение: $$X = 6$$, однако для отыскания его с помощью Maxima придётся прибегнуть к замене исходного уравнения его следствием:

(%i3)	solve([(x-2)=(x-4)^2],[x]);

$$(\%o3)\quad [x=6,x=3]$$

Решения для дальнейшего использования можно извлечь из списка функцией $$ev$$:

(%i1)	sol:solve([x-2=(x-4)^2],[x]);

$$(\%o1)\quad [x=6,x=3]$$

(%i2)	ev(x,sol[1]);

$$(\%o2)\quad 6$$

(%i3)	ev(x,sol[2]);

$$(\%o3)\quad 3$$

Ещё два примера решения алгебраических уравнений:

(%i1)	eq:7*(x+1/x)-2*(x^2+1/x^2)=9;

$$(\%o1)\quad 7\,\left( x+\frac{1}{x}\right) -2\,\left( {x}^{2}+\frac{1}{{x}^{2}}\right) =9$$

(%i2)	sol:solve([eq],[x]);

$$(\%o2)\quad [x=2,x=\frac{1}{2},x=-\frac{\sqrt{3}\,i-1}{2},x=\frac{\sqrt{3}\,i+1}{2}]$$

(%i3)	x1:ev(x,sol[1]); x2:ev(x,sol[2]);
/*комплексные корни не рассматриваем*/

$$(\%o4)\quad 2\frac{1}{2}$$

Уравнения с радикалами перед решением в Maxima приходится преобразовывать к степенной форме (для выделения левой и правой части выражения используют функции $$lhs$$ и $$rhs$$ соответственно):

(%i1)	eq:sqrt(x+1)+sqrt(4*x+13)=sqrt(3*x+12);

$$(\%o1)\quad \sqrt{4\,x+13}+\sqrt{x+1}=\sqrt{3\,x+12}$$

(%i2)	eq1:lhs(eq)^2=rhs(eq)^2;

$$(\%o2)\quad {\left( \sqrt{4\,x+13}+\sqrt{x+1}\right) }^{2}=3\,x+12$$

(%i3)	solve([eq1],[x]);

$$(\%o3)\quad [x=-\sqrt{x+1}\,\sqrt{4\,x+13}-1]$$

(%i4)	eq2:x+1=rhs(%[1])+1;

$$(\%o4)\quad x+1=-\sqrt{x+1}\,\sqrt{4\,x+13}$$

(%i5)	eq3:lhs(eq2)^2=rhs(eq2)^2;

$$(\%o5)\quad {\left( x+1\right) }^{2}=\left( x+1\right) \,\left( 4\,x+13\right)$$

Последняя команда позволила получить степенное уравнение, разрешимое аналитически в Maxima (для этого потребовалось дважды возвести в квадрат исходное уравнение).

(%i6)	solve([eq3],[x]);

$$(\%o6)\quad [x=-4,x=-1]$$

Проверку решения выполняем при помощи функции $$ev$$.

Решение $$x = -4$$ не удовлетворяет исходному уравнению.

(%i7)	ev(eq,%[1]);

$$(\%o7)\quad 2\,\sqrt{3}\,i=0$$

Решение $$x = -1$$ превращает исходное уравнение в верное равенство:

(%i8)	ev(eq,%o6[2]);

$$(\%o8)\quad 3=3$$

Рассмотрим ещё один пример, иллюстрирующий замену и подстановку при решении алгебраических уравнений:

(%i1)	eq:sqrt(x+3-4*sqrt(x-1))+sqrt(x+8-6*sqrt(x-1))=1;

Исходное уравнение:

$$(\%o1)\quad \sqrt{x-4\,\sqrt{x-1}+3}+\sqrt{x-6\,\sqrt{x-1}+8}=1$$

Выполним замену $$\sqrt{x+1}=z,\ z=x^2 + 1$$

(%i2)	eq1:subst(z,sqrt(x-1),eq);

$$(\%o2)\quad \sqrt{-4\,z+x+3}+\sqrt{-6\,z+x+8}=1$$

(%i3)	eq2:subst(z^2+1,x,eq1);

$$(\%o3)\quad \sqrt{{z}^{2}-4\,z+4}+\sqrt{{z}^{2}-6\,z+9}=1$$

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

(%i4)	radcan(%);

$$(\%o4)\quad 2\,z-5=1$$

(%i5)	solve([%],z);

$$(\%o5)\quad [z=3]$$

(%i6)	solve([sqrt(x-1)=3],[x]);

$$(\%o6)\quad [x=10]$$

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

(%i7) ev(eq,%[1]);

$$(\%o7)\quad 1=1$$

Значительная часть тригонометрических уравнений школьного курса также разрешимы в Maxima, но непосредственное решение удаётся получить далеко не всегда.

Примеры:

(%i1)	solve([sin(%pi/6-x)=sqrt(3)/2],[x]);
solve: using arc-trig functions to get a solution.
Some solutions will be lost.

$$(\%o1)\quad [x=-\frac{\pi }{6}]$$

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

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

7.2 Некоторые физические задачи

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

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

7.2.1 Вычисление средней квадратичной скорости молекул

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

Рассмотрим возможность вычисления среднеквадратичной скорости молекул для различных газов. Используемая формула: $$v = \sqrt{\frac{3RT}{M}}$$ , где $$R = 8.314 Дж/(моль \cdot К), T$$ — абсолютная температура, $$M$$ — молярная масса.

Вычислим среднеквадратичную скорость молекул $$CO_2$$ ($$M$$ = 0.044 кг/моль) при температуре 273 К:

(%i1)	v:sqrt(3*R*T/M);

$$(\%o1)\quad \sqrt{3}\,\sqrt{\frac{R\,T}{M}}$$

(%i2)	vCO2:float(v),M=0.044,T=273,R=8.314;

$$(\%o2)\quad 393.3875604633079$$

Расчёт для нескольких различных газов несложно провести, варьируя молярную массу:

(%i3)	vVozd:float(v),M=0.029,T=273,R=8.314;

$$(\%o3)\quad 484.5604478145187$$

(%i4)	vH2:float(v),M=0.002,T=273,R=8.314;

$$(\%o4)\quad 1845.151213315592$$

7.2.2 Распределение Максвелла

Аналогично предыдущему расчёту создадим выражение, описывающее распределение Максвелла (см. блок команд Maxima ниже).

(%i1)	fun:4*%pi*(M/2/%pi/R/T)^(3/2)*exp(-M*v^2/2/R/T)*v^2;

$$(\%o1)\quad \frac{2\,{v}^{2}\,{\left( \frac{M}{R\,T}\right) }^{\frac{3}{2}}\,{e}^{-\frac{{v}^{2}\,M}{2\,R\,T}}}{\sqrt{2}\,\sqrt{\pi }}$$

Для анализа полученных выражений в формулу распределения Максвелла подставляем только температуру:

(%i2)	fun70:fun,T=70;

$$(\%o2)\quad \frac{{v}^{2}\,{\left( \frac{M}{R}\right) }^{\frac{3}{2}}\,{e}^{-\frac{{v}^{2}\,M}{140\,R}}}{35\,\sqrt{2}\,\sqrt{70}\,\sqrt{\pi }}$$

(%i3)	fun150:fun,T=150;

$$(\%o3)\quad \frac{{v}^{2}\,{\left( \frac{M}{R}\right) }^{\frac{3}{2}}\,{e}^{-\frac{{v}^{2}\,M}{300\,R}}}{375\,\sqrt{2}\,\sqrt{6}\,\sqrt{\pi }}$$

(%i4)	fun300:fun,T=300;

$$(\%o4)\quad \frac{{v}^{2}\,{\left( \frac{M}{R}\right) }^{\frac{3}{2}}\,{e}^{-\frac{{v}^{2}\,M}{600\,R}}}{1500\,\sqrt{2}\,\sqrt{3}\,\sqrt{\pi }}$$ (рис 7.1) Распределение Максвелла по скоростям молекул воздуха для различных температур

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

(%i5)	plot2d([fun70,fun150,fun300],[v,0,1000]),M=0.029,R=8.314;

Можно изучить влияние температуры на форму кривой, а так- же на положение максимума функции распределения. С помощью интегрирования $$f(v)$$ можно посчитать долю молекул, обладающих скоростями в каком-либо интервале, а также определить среднюю и среднюю квадратичную скорости молекул.

Пример:

(%i6)	integrate(v*v*fun,v,0,inf);
Is M\,R\,T positive, negative, or zero?
p;

$$(\%o6)\quad \frac{3\,R\,T}{M}$$

Таким образом, $$\langle v^2 \rangle = \frac{3RT}{M}$$ , откуда среднеквадратичная скорость молекул газа $$\langle \sqrt{v^2} \rangle = v_{sq}= \sqrt{\frac{3RT}{M}}$$ .

7.2.3 Броуновское движение

Наличие генератора случайных чисел дает возможность моделировать движение броуновской частицы.

(рис 7.2) Траектория броуновского движения модельной частицы

Эйнштейн первый рассчитал параметры броуновского движения, показав, что нерегулярное перемещение частиц, взвешенных в жидкости, вызвано случайными ударами соседних молекул, совершающих тепловое движение. В соответствии с теорией Смолуховского- Эйнштейна, среднее значение квадрата смещения броуновской частицы ($$s^2$$) за время $$t$$ прямо пропорционально температуре $$T$$ и обратно пропорционально вязкости жидкости $$h$$, размеру частицы $$r$$ и постоянной Авогадро $$N_A:\ s^2 = \frac {2 R T t}{6 \pi h r N_A },$$, где $$R$$ — газовая постоянная.

Броуновские частицы имеют размер порядка 0,1–1 мкм, т.е. от одной тысячной до одной десятитысячной доли миллиметра.

Построим несколько упрощённую модель броуновского движения, предполагая, что смещение частицы по каждой из координат — нормально распределённая случайная величина с нулевым математическим ожиданием. Для генерации случайных чисел используем пакет distrib, включающий необходимые функции (использован генератор $$random_normal$$).

(%i1)	load("distrib")$
	x:0$ y:0$ xy:[[0,0]]$ m:0$ s:1$
	Nmax:500$ for i:1 thru Nmax do (x:x+random_normal(m,s),
	y:y+random_normal(m,s), xy:append(xy,[[x,y]]))$
	plot2d([discrete,xy]);

Результат построения графика приведен на рис. 7.2.

7.3 Пример построения статистической модели

Рассмотрим построение задачи с практическим содержанием.

В табл. 7.1 приведены данные (взяты из статьи В. Ф. Очкова: http://twt.mpei.ac.ru/ochkov/) о зависимости цены подержанного автомобиля от его пробега и "возраста" (времени использования). В статье-первоисточнике задача исследования этой зависимости решалась средствами MathCad.

Рассмотрим её решение средствами Maxima.

Стоимость подержанного автомобиля в зависимости от его возраста и пробега

Возраст (лет)	Пробег (миль)	Цена ($)
11.5	88000	1195
10.5	82000	1295
12.5	97000	800
8.5	51000	2295
9.5	79000	1995
13.5	120000	495
3.5	39000	4995
6.5	52000	2695
4.5	39000	3995
12.5	92000	795
7.5	41000	3495
10.5	77000	1595
12.5	83000	895
4.5	38000	3990
13.5	92000	795
13.5	103000	750
10.5	65000	1495
10.5	70000	1495
10.5	80000	1495
6.5	57000	2695
11.5	101000	895
10.5	78000	1295
9.5	84000	1995
4.5	46000	3675
11.5	108000	975
13.5	124000	850
6.5	56000	3495
9.5	67000	2495
6.5	43000	3400
11.5	78000	1295

В дальнейшем предполагается, что исходные данные для решения подготовлены в виде файла cars.txt. Для считывания используем пакет numericalio. В памяти данные представляются матрицей, а для построения отдельных графиков — списками (переменные $$age, mile, price$$ — см. ниже).

(%i1)	load("draw")$
(%i2)	load("numericalio")$
(%i3)	data:read_matrix("cars1.txt")$
(%i4)	age:makelist(data[k,1], k, 1, 30)$
(%i5)	mile:makelist(data[k,2], k, 1, 30)$
(%i6)	price:makelist(data[k,3], k, 1, 30)$

Простейшую линейную регрессию можно построить, используя функцию $$simple_linear_regression$$ (пакет stats). Построим зависимость цены автомобиля от его стоимости и пробега:

(%i21)	xy:makelist([age[k],price[k]], k, 1, 30)$
(%i22)	simple_linear_regression(xy);

$$(\%o22)\ \begin{pmatrix}SIMPLE\ LINEAR\ REGRESSION\cr model=5757.594446543255-392.7181715149224\,x\cr correlation=-.9688177942467208\cr v\_estimation=95364.34912839333\cr b\_conf\_int=[-431.5987157329751,-353.8376272968697]\cr hypotheses=H0: b = 0 ,H1: b \symbol{`\#} 0\cr statistic=20.69021212080514\cr distribution=[student\_t,28]\cr p\_value=0.0\end{pmatrix}$$

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

(%i26)	xy:makelist([mile[k],log(price[k])], k, 1, 30)$
(%i27)	simple_linear_regression(xy);

$$(\%o27)\ \begin{pmatrix}SIMPLE\ LINEAR\ REGRESSION\cr model=9.174960600286802-2.3747715120748164\,{10}^{-5}\,x\cr correlation=-.9301125564244438\cr v\_estimation=.05467789749118319\cr b\_conf\_int=[-2.7377780810631264\,{10}^{-5},-2.0117649430865062\,{10}^{-5}]\cr hypotheses=H0: b = 0 ,H1: b \symbol{`\#} 0\cr statistic=13.40058098403749\cr distribution=[student\_t,28]\cr p\_value=5.928590951498336\,{10}^{-14}\end{pmatrix}$$

Полученные зависимости представлены в виде выражений Maxima:

(%i28)	fun1:5757.6-392.7*x$(%i29) exp(9.175);

$$9652.768071616591\leqno{(\%o29) }$$

(%i30)	fun2:9653*exp(-2.375*10^(-5)*x);

$$9653\,{e}^{-2.3750000000000001\,{10}^{-5}\,x}\leqno{(\%o30) }$$

Проиллюстрируем полученные результаты графически ( рис. 7.3 и рис. 7.4):

(рис 7.3) Зависимость цены подержанного автомобиля от его возраста (рис 7.4) Зависимость цены подержанного автомобиля от его пробега

(%i34)	draw2d(terminal=eps,key="Table",xlabel="Age",ylabel="Price",
	point_size = 3,point_type=3,points(age,price),
	key="price=f(age)",explicit(fun1,x,0,15));
(%i41)	draw2d(terminal=eps,key="Table",xlabel="Mile",ylabel="Price",
	point_size = 3,point_type=3,points(mile,price),
	key="price=f(mile)",explicit(fun2,x,0,125000));

Для построения модели в виде зависимости цены автомобиля от пробега и возраста одновременно целесообразно использовать более сложную функцию $$lsquares_estimates$$ (пакет lsquares). Искомая модель была представлена уравнением:

$$Price = a + b*Age + c*Mile + d*Mile^2$$

Необходимые команды Maxima:

(%i5)	lsquares_estimates(data,[x,y,z],z=a+b*x+c*y+d*y^2,[a,b,c,d]);

$$(\%o5)\ [[a=\frac{36712000090549571}{5117101479342},b=-\frac{80056614985946}{284283415519},\\ \\ c=-\frac{194393701258481}{3411400986228000},d=\frac{2937180994967}{10234202958684000000}]]$$

(%i6)	float(%);

$$(\%o6)\ [[a=7174.37405506961,b=-281.6084604857839,\\ c=-0.056983539033745,d=2.8699655525931528\,{10}^{-7}]]$$

Следует отметить, что сильно нелинейные задачи решаются при помощи $$lsquares_estimate$$ медленно, поэтому результаты построения модели сильно зависят от обоснованности постановки задачи оценивания. Графическая иллюстрация представлена на рис. 7.5.

(рис 7.5) Иллюстрация зависимости отклика (цены подержанного автомобиля) от двух независимых факторов (возраста и пробега автомобиля)