Решение векторных уравнений в Mathematica

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

A + Vt = Pt

где я знаю A, V и величину P, и мне нужно найти t и направление P. (В принципе, учитывая два луча A и B, где я знаю все об A, кроме начала координат и величины B, выяснить, каким должно быть направление B, чтобы он пересекал A .)

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

Я просмотрел пакет VectorAnalysis и не нашел там ничего подходящего; тем временем пакет Linear Algebra, похоже, имеет решатель только для линейных систем (а это не так, поскольку я не знаю t или P, просто | П|).

Я попытался сделать простодушную вещь: разложить векторы на их компоненты (притвориться, что они трехмерные) и решить их, как если бы я пытался приравнять две параметрические функции,

Solve[ 
      { Function[t, {Bx + Vx*t, By + Vy*t, Bz + Vz*t}][t] == 
          Function[t, {Px*t, Py*t, Pz*t}][t],
        Px^2 + Py^2 + Pz^2 == Q^2 } , 
      { t, Px, Py, Pz } 
     ]

но "решение", которое выплевывается, представляет собой огромную кашу из коэффициентов и заторов. Это также заставляет меня расширять каждое из измерений, которыми я питаю его.

Я хочу красивое символическое решение с точки зрения точечных произведений, перекрестных произведений и норм:

альтернативный текст

Но я не понимаю, как сказать Solve, что некоторые коэффициенты являются векторами, а не скалярами.

Это возможно? Может ли Mathematica дать мне символические решения для векторов? Или я должен просто придерживаться технологии No.2 Pencil?

(Просто для ясности: меня не интересует решение конкретного уравнения вверху — я спрашиваю, могу ли я использовать Mathematica для решения подобных задач вычислительной геометрии, как правило, без необходимости выражать все в виде явной матрицы {Ax, Ay, Az} и др.)


person Crashworks    schedule 21.01.2010    source источник
comment
Возможно, стоит потратить время и задать этот вопрос ребятам из Mathematica. Они, вероятно, знают свое собственное программное обеспечение лучше, чем мы.   -  person Robert Harvey    schedule 21.01.2010
comment
Я попытаюсь разместить сообщение на форумах MathGroup ( forums.wolfram.com ). Не похоже, что простое письмо по адресу [email protected] даст полезные результаты.   -  person Crashworks    schedule 21.01.2010
comment
Удивительно, однако, как часто отношение сигнал/шум (да и вообще вся система) здесь лучше, чем в старомодном списке рассылки/форуме/группе новостей.   -  person Will Robertson    schedule 22.01.2010


Ответы (3)



С Mathematica 7.0.1.0

Clear[A, V, P];
A = {1, 2, 3};
V = {4, 5, 6};
P = {P1, P2, P3};
Solve[A + V t == P, P]

выходы:

{{P1 -> 1 + 4 t, P2 -> 2 + 5 t, P3 -> 3 (1 + 2 t)}}

Ввод P = {P1, P2, P3} может раздражать, если массив или матрица велики.

Clear[A, V, PP, P];
A = {1, 2, 3};
V = {4, 5, 6};
PP = Array[P, 3];
Solve[A + V t == PP, PP]

выходы:

{{P[1] -> 1 + 4 t, P[2] -> 2 + 5 t, P[3] -> 3 (1 + 2 t)}}

Внутренний продукт матричного вектора:

Clear[A, xx, bb];
A = {{1, 5}, {6, 7}};
xx = Array[x, 2];
bb = Array[b, 2];
Solve[A.xx == bb, xx]

выходы:

{{x[1] -> 1/23 (-7 b[1] + 5 b[2]), x[2] -> 1/23 (6 b[1] - b[2])}}

Умножение матриц:

Clear[A, BB, d];
A = {{1, 5}, {6, 7}};
BB = Array[B, {2, 2}];
d = {{6, 7}, {8, 9}};
Solve[A.BB == d]

выходы:

{{B[1, 1] -> -(2/23), B[2, 1] -> 28/23, B[1, 2] -> -(4/23), B[2, 2] -> 33/23}}

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

Однако я не думаю, что перекрестное произведение делает это. Вот как вы используете пакет Notation для его создания. «X» станет нашей инфиксной формой Креста. Я предлагаю скопировать пример из учебника Notation, Symbolize и InfixNotation. Также используйте палитру нотации, которая помогает абстрагироваться от синтаксиса Box.

Clear[X]
Needs["Notation`"]
Notation[x_ X y_\[DoubleLongLeftRightArrow]Cross[x_, y_]]
Notation[NotationTemplateTag[
  RowBox[{x_,  , X,  , y_,  }]] \[DoubleLongLeftRightArrow] 
  NotationTemplateTag[RowBox[{ , 
RowBox[{Cross, [, 
RowBox[{x_, ,, y_}], ]}]}]]]
{a, b, c} X {x, y, z}

выходы:

{-c y + b z, c x - a z, -b x + a y}

Вышеприведенное выглядит ужасно, но при использовании палитры обозначений это выглядит так:

Clear[X]
Needs["Notation`"]
Notation[x_ X y_\[DoubleLongLeftRightArrow]Cross[x_, y_]]
{a, b, c} X {x, y, z}

Я столкнулся с некоторыми причудами, используя пакет обозначений в прошлых версиях mathematica, поэтому будьте осторожны.

person Davorak    schedule 21.01.2010
comment
В Solve[A + V t == P, P] вы пропустили t умножение P. - person rcollyer; 22.01.2010
comment
Вы вводите векторное произведение с :esc: cross :esc:. - person kennytm; 22.01.2010
comment
@rcollyer Моя ошибка, с добавлением должно работать нормально. @KennyTM Спасибо, что указали на это. - person Davorak; 25.01.2010

У меня нет общего решения для вас ни в коем случае (MathForum может быть лучшим способом), но есть несколько советов, которые я могу вам предложить. Во-первых, более систематически разложить ваши векторы на компоненты. Например, я бы решил уравнение, которое вы написали, следующим образом.

rawSol = With[{coords = {x, y, z}},
  Solve[
    Flatten[
     {A[#] + V[#] t == P[#] t & /@ coords,
     Total[P[#]^2 & /@ coords] == P^2}],
    Flatten[{t, P /@ coords}]]];

Тогда вам будет легче работать с переменной rawSol. Далее, поскольку вы ссылаетесь на векторные компоненты единообразно (всегда в соответствии с шаблоном Mathematica v_[x|y|z]), вы можете определить правила, которые помогут их упростить. Я немного поиграл, прежде чем придумал следующие правила:

vectorRules =
  {forms___ + vec_[x]^2 + vec_[y]^2 + vec_[z]^2 :> forms + vec^2,
   forms___ + c_. v1_[x]*v2_[x] + c_. v1_[y]*v2_[y] + c_. v1_[z]*v2_[z] :>
     forms + c v1\[CenterDot]v2};

Эти правила упростят отношения для векторных норм и скалярных произведений (перекрестные произведения оставлены как вероятное мучительное упражнение для читателя). EDIT: rcollyer указал, что вы можете сделать c необязательным в правиле для скалярных произведений, поэтому вам нужно только два правила для норм и скалярных произведений.

С помощью этих правил я сразу смог упростить решение для t до формы, очень близкой к вашей:

  In[3] := t /. rawSol //. vectorRules // Simplify // InputForm
  Out[3] = {(A \[CenterDot] V - Sqrt[A^2*(P^2 - V^2) + 
                                   (A \[CenterDot] V)^2])/(P^2 - V^2), 
            (A \[CenterDot] V + Sqrt[A^2*(P^2 - V^2) + 
                                   (A \[CenterDot] V)^2])/(P^2 - V^2)}

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

person Pillsy    schedule 21.01.2010
comment
Два правила скалярного произведения можно объединить в одно, заменив c_ в первом правиле на c_. или c:1, где точка указывает системе Mathematica использовать стандартное значение по умолчанию для умножения. - person rcollyer; 22.01.2010

У меня несколько иной подход к этому вопросу. Я сделал несколько определений, которые возвращают этот вывод: vExpand examplesМогут быть указаны шаблоны, известные как векторные величины. при использовании vec[_] шаблоны с оберткой OverVector[] или OverHat[] (символы с вектором или шляпой над ними) по умолчанию считаются векторами.

Определения являются экспериментальными и должны рассматриваться как таковые, но, похоже, они работают хорошо. Я ожидаю добавить к этому со временем.

Вот определения. Необходимо вставить в ячейку блокнота Mathematica и преобразовать в StandardForm, чтобы они отображались правильно.

Unprotect[vExpand,vExpand$,Cross,Plus,Times,CenterDot];

(* vec[pat] determines if pat is a vector quantity.
vec[pat] can be used to define patterns that should be treated as vectors.
Default: Patterns are assumed to be scalar unless otherwise defined *)
vec[_]:=False;

(* Symbols with a vector hat, or vector operations on vectors are assumed to be vectors *)
vec[OverVector[_]]:=True;
vec[OverHat[_]]:=True;

vec[u_?vec+v_?vec]:=True;
vec[u_?vec-v_?vec]:=True;
vec[u_?vec\[Cross]v_?vec]:=True;
vec[u_?VectorQ]:=True;

(* Placeholder for matrix types *)
mat[a_]:=False;

(* Anything not defined as a vector or matrix is a scalar *)
scal[x_]:=!(vec[x]\[Or]mat[x]);
scal[x_?scal+y_?scal]:=True;scal[x_?scal y_?scal]:=True;

(* Scalars times vectors are vectors *)
vec[a_?scal u_?vec]:=True;
mat[a_?scal m_?mat]:=True;

vExpand$[u_?vec\[Cross](v_?vec+w_?vec)]:=vExpand$[u\[Cross]v]+vExpand$[u\[Cross]w];
vExpand$[(u_?vec+v_?vec)\[Cross]w_?vec]:=vExpand$[u\[Cross]w]+vExpand$[v\[Cross]w];
vExpand$[u_?vec\[CenterDot](v_?vec+w_?vec)]:=vExpand$[u\[CenterDot]v]+vExpand$[u\[CenterDot]w];
vExpand$[(u_?vec+v_?vec)\[CenterDot]w_?vec]:=vExpand$[u\[CenterDot]w]+vExpand$[v\[CenterDot]w];

vExpand$[s_?scal (u_?vec\[Cross]v_?vec)]:=Expand[s] vExpand$[u\[Cross]v];
vExpand$[s_?scal (u_?vec\[CenterDot]v_?vec)]:=Expand[s] vExpand$[u\[CenterDot]v];

vExpand$[Plus[x__]]:=vExpand$/@Plus[x];
vExpand$[s_?scal,Plus[x__]]:=Expand[s](vExpand$/@Plus[x]);
vExpand$[Times[x__]]:=vExpand$/@Times[x];

vExpand[e_]:=e//.e:>Expand[vExpand$[e]]

(* Some simplification rules *)
(u_?vec\[Cross]u_?vec):=\!\(\*OverscriptBox["0", "\[RightVector]"]\);
(u_?vec+\!\(\*OverscriptBox["0", "\[RightVector]"]\)):=u;
0v_?vec:=\!\(\*OverscriptBox["0", "\[RightVector]"]\);

\!\(\*OverscriptBox["0", "\[RightVector]"]\)\[CenterDot]v_?vec:=0;
v_?vec\[CenterDot]\!\(\*OverscriptBox["0", "\[RightVector]"]\):=0;

(a_?scal u_?vec)\[Cross]v_?vec :=a u\[Cross]v;u_?vec\[Cross](a_?scal v_?vec ):=a u\[Cross]v;
(a_?scal u_?vec)\[CenterDot]v_?vec :=a u\[CenterDot]v;
u_?vec\[CenterDot](a_?scal v_?vec) :=a u\[CenterDot]v;

(* Stealing behavior from Dot *)
Attributes[CenterDot]=Attributes[Dot];

Protect[vExpand,vExpand$,Cross,Plus,Times,CenterDot];
person Codie CodeMonkey    schedule 07.09.2011