Message Boards Message Boards

GROUPS:

How to calculate with symbolic vectors?

Posted 5 months ago
1023 Views
|
3 Replies
|
2 Total Likes
|

I want to generate a system of differential equations for DSolve. The unknown function p[t] is a 2D-function with parameter t and coordinates {x[t], y[t]}. Think of the movement of a point in a plane.

Thereby I encountered a general problem: As long as p is not represented by a list, I cannot tell WL that p is a vector, not a scalar.

I think there should be a symbol attribute, telling that the symbol is meant to be a nD-vector (or nxm matrice, etc.).

See the following piece of code, where p0={p0x, p0y} is a fixed vector, and p shall be a general 2D-Vector. I want to calculate and use d(p) = Norm(p0 - p) within a larger program:

Block[{p0, p0x = 11, p0y = 12, p, x, y, d},
  p0 = {p0x, p0y};
  Print["p0 = ", p0, ", Attr(p0) = ", Attributes[p0], "\tp = ", p, 
   ", Attr(p) = ", Attributes[p]];
  d = Norm[p0 - p];
  Print["1. d(p) = Norm({p0x,p0y} - p) = ", d];
  p = {x, y}; Print["2. d(p=", p, ") = ", d,
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d];
  Clear[x, y];
  d = Norm[p0 - p];
  Print["3. d(p=", p, ") = Norm({p0x,p0y} - {x,y}) = ", d,
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d];
  ];

p0 = {11,12}, Attr(p0) = {} p = p, Attr(p) = {}

1. d(p) = Norm({p0x,p0y} - p) = Sqrt[Abs[11-p]^2+Abs[12-p]^2]

2. d(p={x,y}) = {Sqrt[Abs[11-x]^2+Abs[12-x]^2],Sqrt[Abs[11-y]^2+Abs[12-y]^2]}   d(p={21,22}) = {Sqrt[181],Sqrt[221]}

3. d(p={x,y}) = Norm({p0x,p0y} - {x,y}) = Sqrt[Abs[11-x]^2+Abs[12-y]^2] d(p={21,22}) = 10 Sqrt[2]

In the first step, I calculate Norm[p0 - p], which gives the correct result for a scalar p, but a wrong result for a vector p. If in the second step I replace p by {x,y} within the result of the first step, I get the wrong result of course.
In the third step, I recalculate Norm[p0 - p], which is correct now.

Of course I can redefine d(p) by a function with set-delayed: d[p_] := Norm[p0 - p] Then everything works fine, if I use d[p] instead of Norm[p0 -p] and replacement.

Block[{p0, p0x = 11, p0y = 12, p, x, y, d},
  p0 = {p0x, p0y};
  Print["p0 = ", p0, ", Attr(p0) = ", Attributes[p0], "\tp = ", p, 
   ", Attr(p) = ", Attributes[p]];
  d[p_] := Norm[p0 - p];
  Print["1. d(p) = Norm({p0x,p0y} - p) = ", d[p]];
  p = {x, y}; Print["2. d(p=", p, ") = Norm({p0x,p0y} - p] = ", d[p],
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d[p]];
  Clear[x, y];
  Print["3. d(p=", p, ") = Norm({p0x,p0y} - {x,y}) = ", d[p],
   {x, y} = {21, 22}; "\td(p=", p, ") = ", d[p]];
  ];

p0 = {11,12}, Attr(p0) = {} p = p, Attr(p) = {}

1. d(p) = Norm({p0x,p0y} - p) = Sqrt[Abs[11-p]^2+Abs[12-p]^2]

2. d(p={x,y}) = Norm({p0x,p0y} - p] = Sqrt[Abs[11-x]^2+Abs[12-y]^2] d(p={21,22}) = 10 Sqrt[2]

3. d(p={x,y}) = Norm({p0x,p0y} - {x,y}) = Sqrt[Abs[11-x]^2+Abs[12-y]^2] d(p={21,22}) = 10 Sqrt[2]

The problem is, that this doesn't help in my original problem of a differential equation for a vector-valued function p[t].

See the following code for a trivial constant-velocity movement in a plane, starting at some point p0:

Block[{p0 = {1, 2}, t, p, v = {2, 1}, eqs, sol},
  eqs = {p'[t] == v, p[0] == p0};
  Print["eqs = ", Column[eqs]];
  sol = DSolve[eqs, p[t], {t, 0, 4}];
  Print["sol = ", sol];
  ];

eqs = (p^\[Prime])[t]=={2,1}
p[0]=={1,2}

DSolve::nolist: List encountered within {(p^\[Prime])[t]=={2,1}}. There should be no lists on either side of the equations.

sol = DSolve[{(p^\[Prime])[t]=={2,1},p[0]=={1,2}},p[t],{t,0,4}]

Even though by the boundary condition p[0]==p0 it is clear that p is a vector, WL doesn't recognize p[t] and therefore p'[t] to be vector-valued functions of a scalar t.

Since I cannot tell WL that p[t] is vector-valued, I am forced to switch to coordinate-based differential equations:

Block[{p0 = {1, 2}, t, p, x, y, v = {2, 1}, eqs, sol},
  p[t_] = {x[t], y[t]};
  eqs = {p'[t] == v, p[0] == p0};
  Print["eqs = ", Column[eqs]];
  sol = DSolve[eqs, p[t], {t, 0, 4}];
  Print["sol = ", sol];
  ];

eqs = {(x^\[Prime])[t],(y^\[Prime])[t]}=={2,1}
{x[0],y[0]}=={1,2}

sol = {{x[t]->1+2 t,y[t]->2+t}}

This seems to be no big issue but imagine p'[t] as a complicated expression in p[t] and t. You have to switch to NDSolve then and you will get a list of scalar InterpolatingFunction(s) as the solution for {x[t], y[t]}. It would be much better to get one vector-valued InterpolationFunction as the solution for p[t].

3 Replies
Posted 4 months ago

Hello Werner,

my prefered way for doing computation with symbolic vectors is the following: I define the vectors using Indexed and Array:

ndim = 2;
p0V = Array[Indexed[p0, #] &, ndim]; 

pV = Array[Indexed[p, #][t] &, ndim]; 
vV = Array[Indexed[v, #] &, ndim];

I use the convention to append "V" to the name of the vector in order to avoid any confusion between the vector and its components. Computation of the Norm is straightforward:

Norm[p0V - pV]

Next the differential equation:

eqs = { D[pV, t] == vV, (pV /. t -> 0) == p0V} 

sol = Flatten[DSolve[eqs, pV, t]]

If you want to have a more comprehensive printout the indices can be replaced

sol /. Indexed[u_, n_] -> Indexed[u, {"x", "y"}[[n]]]

Using ndim and Array you can easily switch from say 2D to 3D models, of course you have to adjust the specification of the boudary condition for the DGL and the number of labels for renaming

Regards, Michael

Posted 4 months ago

Hello Michael,

That's a nice idea to declare vectors as general indexed arrays through:

pV=Array[Indexed[p,#][t]&,ndim]

The advantage is when you have a lot of vectors, a single change in ndim brings them all to a new length. The disadvantage is, it is not easy to read and understand. The definition:

pV[t_]={px[t],py[t]} (* or with subscripted x, y *)

is much more readable.

But more import those indexed arrays don't change anything of the fact that we are switching from vector-equations to coordinate equations. I.e:

DSolve[{pV'[t]==vV, ...}, pV, t]

actually solves the scalar system

DSolve[{{px'[t], py'[t]}=={vx, vy}, ...}, {px[t], py[t]}, t]

and not the original vector-DGL.

That is why I proposed to give symbols a new attribute "2DVector" or "Dimensions" or the like.

BTW: The replacement "Indexed[u,{"x","y"}[[n]]]" works, but gives a strange error message "... Part: The expression n cannot be used as a part specification". I don't understand, why.

It seems to me that at the very start, the expression ndim = 2 defeats the intent of the question as I understand it, namely, to have an arbitrary, symbolic dimension n and not some particular dimension.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract