Maple 9.5/10 в математике, физике и образовании - Владимир Дьяконов
Шрифт:
Интервал:
Закладка:
В пакете векторных операций определен ряд типовых операций с кривыми. Ниже представлено задание эллиптической кривой и вычисление в аналитической форме нормали и радиуса кривизны (файл vopcurves):
> SetCoordinates(cartesian);
cartesian> assume(t::real);
> ell := <2*cos(t),sin(t)>;
ell := 2 cos(t)ex + sin(t)ey> nv := simplify(PrincipalNormal(ell,t));
> len := simplify(LinearAlgebra:-Norm(nv, 2));
> r := simplify(RadiusOfCurvature(ell));
Теперь можно представить саму кривую (эллипс) и ее эволюту (рис. 4.39):
> ev := simplify(ell + r * nv / len);
> plot([[ell[1], ell[2], t=0..2*Pi], [ev[1], ev[2], t=0..2*Pi]]);
Рис. 4.39. Графики кривой — эллипса и ее эволюты
Нетрудно заметить, что для эллипса эволюта представляет собой удлиненную астроиду.
Для вычисления кривизны кривой С используется функция Curvature(C, t) в которой параметр t может и отсутствовать:
> Curvature(<cos(t),t,sin(t)>, t);
> с := Curvature(t -> <t,t^2,t^4>):
> simplify(c(t)) assuming t::real;
> SetCoordinates('polar');
polar> Curvature(<exp(-t^2), t>):
> simplify(%) assuming t::real;
4.11.5. Интегрирование в пакете VectorCalculus
В аспекте практических приложений векторного анализа и теории поля особый интерес представляют приложения интегрирования пакете VectorCalculus. Так, видоизмененная функция int(f, dom) задает вычисление интеграла от функции f по области dom, например (файл vecint):
> restart:with(VectorCalculus):
> int(х^2+у^2, [x,y] = Circle(<0,1>, r));
> int(sin(х)*cos(у)*tan(z), [x,y,z] = Parallelepiped(0..Pi, 0..Pi/3, 0..Pi/4));
½√3 ln(2)Функция PathInt(f, dom) вычисляет интеграл пути для функции f с Rn до R:
> PathInt(х^2, [х,y] = Line(<0,0>, <1,2>));
> PathInt(х^2+y^2, [х,y] = Circle(<0,0>, 3/2));
> PathInt(1, [х,y] = Ellipse(х^2+y^2/2-1));
Другая функция LineInt(F, dom), где F — вектор или процедура задания векторного поля, dom — параметр, характеризующий направление интегрирования, задает вычисление линейного интеграла в пространстве Rn:
> SetCoordinates(cartesian[х,y]);
cartesianx, у> LineInt(VectorField(<х,y>), Line(<0,1>, <2,-5>));
14> LineInt(VectorField(<y,-х>), Circlet<0,0>, r));
-2 r² π> LineInt(VectorField(<y,-х>), Ellipse(х^2/4+y^2/9-1));
-12π> LineInt(VectorField(<y,-х>), Arc(Ellipse(х^2/4+у^2/9-1), 0, Pi/2));
-3πФункция ArcLength(C,dom) задает вычисление длины дуги С по известному интегральному выражению для нее:
> ArcLength(<r*cos(t),r*sin(t)>, t=0..Pi) assuming r>0;
πr> ArcLength(t -> <t,t^2>, 0..2);
√17-¼ln(-4+√17)> evalf(%);
4.646783762Рекомендуется просмотреть различные варианты задания области интегрирования dom в справке по этому пакету.
4.11.6. Задание матриц специального типа
Пакет VectorCalculus позволяет для заданной функции f задавать несколько матриц специального вида, которые часто используются при решении задач теории поля:
Hessian(f, t) — создание матрицы гессиана;
Jacobian(f, v, det) — создание матрицы якобиана;
Wronskian(f, t) — создание матрицы вронскиана.
Примеры задания таких матриц приведены ниже (файл vecmatrix):
> Hessian(ехр(х*y), [х,y]);
> Hessian(а/(х^2+y^2+z^2), [х, y, z]);
> Н := unapply(%, [a,x,y,z]):
> Н(1/2, 0.3, 0.7, 0.1);
> Jacobian([r*cos(t), r*sin(t)], [r,t]);
> Jacobian([r*cos(t), r*sin(t)], [r,t], 'determinant');
> Wronskian([exp(t),ln(t),sin(t)], t);
> Wronskian([t, t^2, t^3], t)
4.11.7. Функции теории поля
К основным функциям теории поля относятся:
Curl(F) — вычисляет вихрь векторного поля в R³;
Divergence(F) — вычисляет дивергенцию векторного поля;
Flux(f, dom) — вычисляет поток векторного поля в R³;
Gradient(f, с) — вычисляет градиент функции f в пространстве от Rn до R;
Del(f, с) и Nabla(f, с) — векторные дифференциальные операторы;
Laplacian(f, с) или Laplacian(F) — вычисляет лапласиан функции f или векторного определения (процедуры) F;
ScalarPotential(v) — вычисляет скалярный потенциал векторного поля;
Torsion(C, t) — вычисляет торсион в R³;
VectorPotential(v) — вычисляет векторный потенциал в R³;
Довольно громоздкие определения этих функций, основанные на использовании криволинейных и поверхностных интегралов, имеются в учебной литературе. Не приводя их, ограничимся приведенными ниже примерами применение указанных выше функций (файл vecft):
> restart:with(VectorCalculus): SetCoordinates('cartesian'[x,y,z]);
cartesianx, у, z> F := VectorField( <-y,x,0> );
F:=-yēx +хēу> Curl(F);
2ēz> Del &x F;
2ēz> Nabla &x F;
2ē> CrossProduct(Del, F);
2ēz> F := VectorField(<х^2,y^2,z^2>);
F:=-x²ēх + y²ēу + z²ēz> Divergence(F);
2х + 2у + 2z> Flux(VectorField(<x,y,z>, cartesian[x,y,z]), Sphere(<0,0,0>, r));
4r³ π> Gradient(х^3/3+у^2, [x,y]);
x²ēx + 2yēу 0ēх> Del(х^2+у^2+z^2);
2xēx + 2уēу + 2zēz> Nabla(х^2+у^2+z^2);
2xēx + 2уēу + 2zēz> Del . %;
6> Laplacian(х^2+у^2+z^2, [x,y,z]);
6> Laplacian(f(r,theta,z));
> SetCoordinates('cylindrical' [r, theta, z])
cylindricalr, θ, z> Laplacian(f(r, theta, z));
> SetCoordinates('cartesian'[x,y,z]);
cartesianx, y, z> v := VectorField(<x,y,-z>);
v := xēx + уēу - zēz> ScalarPotential(v);
> v := VectorField(<-y,0,z>);
v := -yēx + zēz> ScalarPotential(v); den := х^2 + y^2 + z^2;
den := x² + y² + z²> ScalarPotential((x,y,z) -> <x,y,z>/den);
(x,y,z)→½ ln(x² + y² + z²)> SetCoordinates('spherical'[r,phi,theta]);
sphericalr, φ, θ> v := VectorField(<r,0,0>);
v:= r ēг> ScalarPotential(v);
> restart:with(VectorCalculus): simplify( Torsion(<t,t^2,t^3>)) assuming t::real;
> Torsion(t -> <2*t,sin(t),cos(t)>);
> SetCoordinates('cartesian'[x,y,z]); v := VectorField(<y,-x,0>);
cartesianx, y, z v:= уēx - хēу> VectorPotential(v);
-xzēx - yzēу> SetCoordinates('cylindrical'[r,theta,z]);
cylindricalr, θ, z> v := VectorField(<r,0,-2*z>);
v:= rēr -2zēz> VectorPotential(v);
(-r sin(θ)² z - r cos(θ)² z) ēθ> simplify(Curl(%));
rēr - 2zēzОбратите внимание на то, что для гарантии правильного выполнения этих команд и отсутствия «зависания» компьютера может потребоваться команда restart и перезагрузка пакета VectorCalculus.
4.11.8. Приближение площади сложной поверхности суммами Римана
Одним из важнейших приложений пакета VectorCalculus является вычисление длин дуг и площадей сложных поверхностей на основе применения линейных и поверхностных интегралов. Иногда это встречает большие трудности и требует специальных подходов. Примером может служить поверхность, заданная рис. 4.40. Эта поверхность построена с имитацией ее освещения от внешнего источника света.
Рис. 4.40 Сложная поверхность с эффектами ее освещения внешним источником света
Применим обычную процедуру вычисления площади поверхности. Для этого вычислим для нее матрицу якобиана и удалим из нее столбец с нулевыми элементами (файл vecrim):
> J := Jacobian(f, [х, у, z]);
> J := DeleteColumn(J, [3]);
Тогда площадь поверхности вычисляется следующим образом:
> Int(Int(dA, x=0..2*Pi), y=0..2*Pi);
К сожалению, этот двойной интеграл Maple не вычисляет из-за сложности подынтегрального выражения, график которого представлен на рис. 4.41.
Рис. 4.41. График подынтегрального выражения
Для приближенного вычисления площади можно разбить поверхность на достаточное число сегментов и использовать замену интегралов суммами Римана. Оценка нижней и верхней сумм Римана для четверти поверхности (ее одного квадранта) представлена ниже:
