A motion profile problem or Mathcad sucks. Will some other CAS solve this?
pnachtwey
an initial PVA ( position, velocity and acceleration ) to a final PVA
where the final V=0 and A=0. Sound easy, I have 6 equations and 6
unknows to solve for. There are three motion segments. The first is
from point 0 to point 1 where the jerk is postive and the target is
accelerating. The second is from point 1 to point 5. ( i skipped
some points ) where the jerk is negative and a5 is at peak
deceleration. The final segment goes from point 5 to point 7. Duing
this segment the jerk is postive again and the deceleration is
decreasing until the target stops at x7.
Here is a link to my Mathcad work sheet.
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Seg137.pdf
It looks like it should be simple. Mathcad does come up with an answer
but after that it doesn't respond to the keyboard or mouse very well
or at all. I keep getting a Not Responding errorr. I can't cut and
paste the resulting RootOf equation so I can't solve for the roots and
do the manual operations. BTW, I hate having to cut and paste
equations or intermediate results when I change a parameter.
Will Mathematica, Maple or other CAS solve this, symbolically,
easily, quickly? Can it produce equations in text form so I can cut
and paste it into my C program? At this time I must retype all the
Mathcad formulas into C. This is tedious and error prone.
While I am at it I would like is a semi automatic mode. Many times I
am smarter than the Mathcad. I would like to be able to do simple
functions like add, subtract, multiply or divide by the same
expression on both side of equation easily. I would like to multiply
the numerator or denominator of a fraction by the same number easily.
Perhaps this mode wouldn't be necesary in a better CAS package.
I would like a CAS that I could use to write a book or manual.
Mathcad does a faily good job in this area but it has the nasty habbit
of display a-b as -b+a all the time and the .rtf formating is
sometimes poor. Mathcad will not generate html files that work on
unix/linux servers. It uses \\ instead of // in the URL and file
names so images don't show up when using a unix server.
BTW, I have been lurking this group for a while now. Vladimir B has
me wondering if any CAS good enough. I am an engineer. I generate a
lot of symbolic equations for C, C++ programs but I don't get into the
far out stuff that you Mathematicians do. I optimize/minimuze, solve
differential equations and use Laplace and z transforms. What would
be cool is too solve some of the simpler differential equations
symbolically instead of using Runge Kutta. Are there any CAS packages
that can at least to that without running into bugs? I can do the
math I need. I just need a tool that works.
Does anything fit my needs? There must be. I just don't want to buy
a bunch of bugs. Mathcad is pitiful compared to the stuff I see on
this group
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20pitiful.pdf
Mathcad13 crashes and hangs all the time up but at least I have not
found any wrong answers, yet.
Has anybody got a Mathcad13 to Mathematica or Maple converter that
works?
I know this is getting to be more like a wish list.
Thanks in advance
Peter Nachtwey
Roman Pearce
> ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Seg137.pdf
> Will Mathematica, Maple or other CAS solve this, symbolically,
> easily, quickly? Can it produce equations in text form so I can cut
> and paste it into my C program? At this time I must retype all the
> Mathcad formulas into C. This is tedious and error prone.
I used the numerators of the equations (ignoring that j != 0).
Maple 11 will compute a lexicographic Groebner basis for the system
where {x0,v0,a0,x7,j} are parameters and {a1,x1,a5,x5,v5,v1} are
variables using the F4 and FGLM algorithms. This is one of those
cases where it pays to implement algorithms in general domains.
Time was 0.37 sec on a 1.4 GHz Athlon.
I hope I copied down the equations correctly:
eqns := {
x0 + v0*(a1-a0)/j + a0/2*(a1-a0)^2/j^2 + j/6*(a1-a0)^3/j^3 - x1,
v0 + a0*(a1-a0)/j + j/2*(a1-a0)^2/j^2 - v1,
x1 + v1*(a1-a5)/j + a1/2*(a1-a5)^2/j^2 - j/6*(a1-a5)^3/j^3 - x5,
v1 + a1*(a1-a5)/j - j/2*(a1-a5)^2/j^2 - v5,
x5 + v5*(-a5/j) + a5/2*(-a5/j)^2 + j/6*(-a5/j)^3 - x7,
v5 + a5*(-a5/j) + j/2*(-a5/j)^2
};
sys := map(numer@expand, eqns);
vars := {x1,v1,a1,x5,v5,a5};
G := Groebner[Basis](sys, 'tord', variables=vars, order=plex):
Here are the polynomials with x1 > x5 > v1 > v5 > a1 > a5.
You will need to solve the first polynomial for a5 in terms of the
parameters, then substitute that into the second polynomial and
solve for a1, etc. I wasn't able to factor the system at all,
so a general "symbolic" solution is going to be messy.
[
144*x0*j^4*x7-72*x0^2*j^4-72*x7^2*j^4+144*x0*j^3*v0*a0-144*v0*a0*j^3*x7
-72*v0^3*j^3+36*v0^2*j^2*a0^2+48*x7*j^2*a0^3-48*x0*j^2*a0^3-6*v0*j*a0^4
+a0^6+(72*v0*j*a0^2-72*v0^2*j^2-18*a0^4)*a5^2+(-144*x7*j^2+144*x0*j^2
+48*a0^3-144*v0*a0*j)*a5^3+(72*v0*j-36*a0^2)*a5^4,
(72*v0^3*j^3+72*x0^2*j^4+72*x7^2*j^4-a0^6+48*x0*j^2*a0^3-144*x0*j^4*x7
-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*x7*j^2*a0^3-144*x0*j^3*v0*a0
+144*v0*a0*j^3*x7)*a1+(72*v0^3*j^3+72*x0^2*j^4+72*x7^2*j^4-a0^6
+48*x0*j^2*a0^3-144*x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*x7*j^2*a0^3
-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7)*a5+(36*a0^2*x0*j^2+72*v0^2*j^2*a0
-72*v0*j^3*x0+12*a0^5-60*v0*j*a0^3+72*v0*j^3*x7-36*a0^2*x7*j^2)*a5^2
+(72*v0*j*a0^2-72*v0^2*j^2-18*a0^4)*a5^3,
2*v5*j-a5^2,
-2*v0*j+a0^2-2*a5^2+4*v1*j,
6*x5*j^2-6*x7*j^2-a5^3,
-1080*x0^3*j^6+5*a0^9-216*x7^3*j^6+2268*v0^2*j^4*a0^2*x7
-1458*v0*j^3*a0^4*x7+1350*x0*j^3*a0^4*v0+1296*x0*j^4*a0^3*x7
-1620*v0^2*j^4*a0^2*x0-3888*x0*j^5*x7*v0*a0+648*x7^2*j^5*v0*a0
+3240*x0^2*j^5*v0*a0-1080*v0^3*j^5*x0+1080*v0^4*j^4*a0-900*v0^3*j^3*a0^3
-216*v0^3*j^5*x7-1080*x0^2*j^4*a0^3+1944*x0^2*j^6*x7-648*x7^2*j^6*x0
-216*x7^2*j^4*a0^3-225*a0^6*x0*j^2-45*a0^7*v0*j+243*a0^6*j^2*x7
+270*v0^2*j^2*a0^5+(864*j^5*v0*x0^2+864*j^5*v0*x7^2-1728*j^5*v0*x0*x7
+864*j^4*v0^4-1728*j^4*v0^2*x0*a0+1728*j^4*v0^2*a0*x7-432*j^4*x7^2*a0^2
-432*j^4*x0^2*a0^2+864*j^4*a0^2*x0*x7-864*j^3*v0^3*a0^2+1440*j^3*v0*x0*a0^3
-1440*j^3*v0*x7*a0^3+288*j^2*v0^2*a0^4-288*j^2*x0*a0^5+288*j^2*a0^5*x7
-48*v0*a0^6*j+6*a0^8)*a5+
(-864*x0*j^4*v0^2+864*x7*j^4*v0^2+864*v0^3*a0*j^3
-864*x7*j^3*v0*a0^2+864*x0*j^3*v0*a0^2-1152*v0^2*a0^3*j^2+216*x7*j^2*a0^4
-216*x0*j^2*a0^4+504*v0*a0^5*j-72*a0^7)*a5^2+(432*x0*j^4*x7-216*x0^2*j^4
-216*x7^2*j^4+432*x0*j^3*v0*a0-432*v0*a0*j^3*x7-1080*v0^3*j^3
-144*x0*j^2*a0^3+1404*v0^2*j^2*a0^2+144*x7*j^2*a0^3-666*v0*j*a0^4
+111*a0^6)*a5^3+(-2592*j^5*v0*x0*a0+2592*j^5*v0*a0*x7+1296*v0^3*j^5
-2592*j^6*x0*x7-648*v0^2*j^4*a0^2+864*j^4*a0^3*x0-864*j^4*a0^3*x7
+108*v0*j^3*a0^4+1296*j^6*x0^2+1296*j^6*x7^2-18*a0^6*j^2)*x1
];
Peter Nachtwey
"Roman Pearce" <rpea...@gmail.com> wrote in message
news:06e1982a-5435-48a8...@q27g2000prf.googlegroups.com...
> pnachtwey wrote:
>> ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Seg137.pdf
>> Will Mathematica, Maple or other CAS solve this, symbolically,
>> easily, quickly? Can it produce equations in text form so I can cut
>> and paste it into my C program? At this time I must retype all the
>> Mathcad formulas into C. This is tedious and error prone.
>
> I used the numerators of the equations (ignoring that j != 0).
a=acceleration and j=jerk or third derivative. The jerk can't be zero is
the rate of change in acceleration. A typical value is about 10g/s. I
normally use 4000 in/s^3 or 100 m/s^3 depending on the unit I am using.
> Maple 11 will compute a lexicographic Groebner basis for the system
> where {x0,v0,a0,x7,j} are parameters and {a1,x1,a5,x5,v5,v1} are
> variables using the F4 and FGLM algorithms.
This is all new to me. I have some studying to do. I did quickly look
these algroithms up on Wikipedia but there isn't enough back ground there
figure out exactly what they are taling about.
> This is one of those
> cases where it pays to implement algorithms in general domains.
> Time was 0.37 sec on a 1.4 GHz Athlon.
It took at least 2 minutes on my 3.4 GHz hyper threaded laptop.
>
> I hope I copied down the equations correctly:
The equation look ok if you are assuning they must equate to 0
>
> eqns := {
> x0 + v0*(a1-a0)/j + a0/2*(a1-a0)^2/j^2 + j/6*(a1-a0)^3/j^3 - x1,
> v0 + a0*(a1-a0)/j + j/2*(a1-a0)^2/j^2 - v1,
> x1 + v1*(a1-a5)/j + a1/2*(a1-a5)^2/j^2 - j/6*(a1-a5)^3/j^3 - x5,
> v1 + a1*(a1-a5)/j - j/2*(a1-a5)^2/j^2 - v5,
> x5 + v5*(-a5/j) + a5/2*(-a5/j)^2 + j/6*(-a5/j)^3 - x7,
> v5 + a5*(-a5/j) + j/2*(-a5/j)^2
> };
> sys := map(numer@expand, eqns);
> vars := {x1,v1,a1,x5,v5,a5};
> G := Groebner[Basis](sys, 'tord', variables=vars, order=plex):
This is interesting.
>
> Here are the polynomials with x1 > x5 > v1 > v5 > a1 > a5.
> You will need to solve the first polynomial for a5 in terms of the
> parameters, then substitute that into the second polynomial and
> solve for a1, etc. I wasn't able to factor the system at all,
> so a general "symbolic" solution is going to be messy.
I will try to enter this into mathcad and verify.
What I find interesting is that Mathcad returns a complicated
RootOf(polynomial of _Z). I must find the roots and the pick the right one
which is usually the positive real one.Maple seems to pick the right root.
I then must substitute this value for the RootOf expression. This is very
manual. Perhaps all this is automatic in Maple because of the algortihms
you mentioned above.
I don't understand why a5 should be so much more complicated than a1. I
would have thought a5 would be simpler than a1 because the final
acceleration and velocity are 0 whereas the initial velocity and
acceleration ( x0 and a0 ) may not be. At least a symbolic solution is
returned.
The solution for a5 is too complicated to be computed on a motion controller
on-the-fly but so is finding the 4th ,5th or 6th root of a polynomial.
Thanks. I will plug it in. I have many similar Mathcad worksheets. I can
see I have been wasting my time solving these kinds of problems with
Mathcad.
Thanks, the information about the algorthms is probably more important than
the single solution. This has raised the high bar. I will study them.
Peter Nachtwey
pnachtwey
If this is the answer to x1 it can't be right. Noticie that a5 is one
of the parameters. x1, v1, a1, x5, v5, a5 must be solved only in
terms of x0, v0, a0, x7 and j. I see a5 in the above equation.
This means that I must compute a5 first, ugly. Computing a5 first
means I must compute x1 first, but x1 is dependent on a5. This can't
work.
Either I don't understand how to interpret the answer or the answer is
wrong. I would expect to get an equation for each of
x1,v1,a1,x5,v5,a5 in terms of x0, v0, a0, x7 and j. I can see
calculating a simplest term of x1,v1,a1,x5,v5,a5 and substituting it
into the other terms but circular substitution will not work.
Here is an example of a short motion profile where Mathcad does work
after manual manipulation to copy and paste the polynomial in the
RootOf results and finding the roots but the answer is right. For
some reason Mathcad doesn't hang on this problem even though there are
more equations and unknowns.
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Seg13567%20%20LT%20ReqVel.pdf
Notice that the solutions are in terms of x0, v0, a0, x7 and j only
and there are no circular definitions.
Can Mathematica or Maple solve this without making me do the manual
operations?
Peter Nachtwey
Roman Pearce
> Either I don't understand how to interpret the answer or the answer is
> wrong. I would expect to get an equation for each of
> x1,v1,a1,x5,v5,a5 in terms of x0, v0, a0, x7 and j.
For what I computed, you get equations for a5, a1, v5, v1, x5, x1,
where each equation involved the parameters {x0,v0,a0,x7,j} and the
previous variables, ie: the equation for v1 may contain {a5,a1,v5} as
well as {x0,v0,a0,x7,j}. This is a triangular form. You solve for
the lowest variable (a5) first, plug the symbolic solution into the
remaining equations, then solve the next lowest variable (a1), and
repeat. Because the system has a finite number of solutions in
{a5,a1,v5,v1,x5,x1}, this procedure will eventually give you
"solutions" for {a5,a1,v5,v1,x5,x1} expressed as roots of polynomials
in {x0,v0,a0,x7,j}.
There are other ways to arrive at an equivalent result. By
introducing a new variable Z = random linear combination of the
variables, and trying to "separate" the solutions in terms of Z. ie:
get {a5,a1,v5,v1,x5,x1} equal to polynomials in Z, where Z =
RootOf(giant mess). This may not be better.
I think it's fair to say that your system has no nice solutions. Does
anyone else want to take a crack at it ?
pnachtwey
> On Apr 1, 10:37 am, pnachtwey <pnacht...@gmail.com> wrote:
>
> > Either I don't understand how to interpret the answer or the answer is
> > wrong. I would expect to get an equation for each of
> > x1,v1,a1,x5,v5,a5 in terms of x0, v0, a0, x7 and j.
>
> For what I computed, you get equations for a5, a1, v5, v1, x5, x1,
> where each equation involved the parameters {x0,v0,a0,x7,j} and the
> previous variables, ie: the equation for v1 may contain {a5,a1,v5} as
> well as {x0,v0,a0,x7,j}. This is a triangular form. You solve for
> the lowest variable (a5) first, plug the symbolic solution into the
> remaining equations, then solve the next lowest variable (a1), and
> repeat. Because the system has a finite number of solutions in
> {a5,a1,v5,v1,x5,x1}, this procedure will eventually give you
> "solutions" for {a5,a1,v5,v1,x5,x1} expressed as roots of polynomials
> in {x0,v0,a0,x7,j}.
I understand but I don't see any of the 6 solutions without a a5 in
it. This can't be right.
>
> There are other ways to arrive at an equivalent result. By
> introducing a new variable Z = random linear combination of the
> variables, and trying to "separate" the solutions in terms of Z. ie:
> get {a5,a1,v5,v1,x5,x1} equal to polynomials in Z, where Z =
> RootOf(giant mess). This may not be better.
This is what Mathcad appears to use. When it doesn't hang up it
appears to get the right answer even though it isn't really a symbolic
one and requires manual manipulation
>
> I think it's fair to say that your system has no nice solutions. Does
> anyone else want to take a crack at it ?
I would appreciate it. It is only 6 equations and 6 unknowns. How
hard can it be? ;)
Even knowing that CAS can't solve this lets me know that I shouldn't
waste my money on any of them and that is worth knowing too.
Peter Nachtwey
Roman Pearce
Y=RootOf((72*v0*j-36*a0^2)*_Z^4+(-144*x7*j^2-144*v0*a0*j
+144*x0*j^2+48*a0^3)*_Z^3+(-18*a0^4+72*v0*j*a0^2-72*v0^
2*j^2)*_Z^2-72*x7^2*j^4+144*x0*j^4*x7-72*x0^2*j^4+144*x0*j^3*v0*a0-72*v0^3*j^3-144*v0*a0*j^3*x7+36*v0^2*j^2*
a0^2-48*x0*j^2*a0^3+48*x7*j^2*a0^3+a0^6-6*v0*j*a0^4);
That is, Y is a root of a degree 4 polynomial with coefficients in {j,
x0, v0, a0, x7}. Then the solutions are:
[x1 =
1/18*(1080*v0^3*j^5*x0-1080*v0^4*j^4*a0+225*a0^6*x0*j^2+648*x7^2*j^6*x0-243*a0^6*x7*j^2-648*x7^2*j^5*v0
*a0+216*x7^2*j^4*a0^3-270*v0^2*j^2*a0^5+216*v0^3*j^5*x7+1728*Y*v0*j^5*x0*x7-1440*Y*v0*j^3*x0*a0^3+1440*Y*v0*j
^3*x7*a0^3-1728*Y*v0^2*j^4*a0*x7+1728*Y*v0^2*j^4*x0*a0-864*Y*a0^2*x0*j^4*x7+1080*x0^2*j^4*a0^3+216*Y^3*x0^2*j
^4+216*Y^3*x7^2*j^4+1080*Y^3*v0^3*j^3-864*Y*v0^4*j^4-3240*x0^2*j^5*v0*a0-1944*x0^2*j^6*x7+666*Y^3*v0*j*a0^4+
864*Y^2*x0*j^4*v0^2-864*Y^2*v0^3*a0*j^3+1152*Y^2*v0^2*a0^3*j^2-216*Y^2*x7*j^2*a0^4-504*Y^2*v0*a0^5*j-864*Y^2*
x7*j^4*v0^2+216*Y^2*x0*j^2*a0^4+1458*v0*j^3*x7*a0^4-2268*v0^2*j^4*x7*a0^2-111*Y^3*a0^6+72*Y^2*a0^7-6*Y*a0^8+
216*x7^3*j^6-5*a0^9+45*v0*j*a0^7+900*v0^3*j^3*a0^3+1620*v0^2*j^4*x0*a0^2-1350*v0*j^3*x0*a0^4+1080*x0^3*j^6-
\
1296*x0*j^4*a0^3*x7+3888*x0*j^5*x7*v0*a0-432*Y^3*x0*j^3*v0*a0+432*Y^3*v0*a0*j^3*x7+864*Y^2*x7*j^3*v0*a0^2-864
*Y^2*x0*j^3*v0*a0^2+864*Y*v0^3*j^3*a0^2-864*Y*v0*j^5*x0^2+432*Y*a0^2*x7^2*j^4-864*Y*v0*j^5*x7^2+48*Y*v0*j*a0^
6+432*Y*a0^2*x0^2*j^4+288*Y*a0^5*x0*j^2-288*Y*a0^5*x7*j^2-288*Y*v0^2*j^2*a0^4-432*Y^3*x0*j^4*x7-1404*Y^3*v0^2
*j^2*a0^2+144*Y^3*x0*j^2*a0^3-144*Y^3*x7*j^2*a0^3)/j^2/
(72*x0^2*j^4+72*x7^2*j^4+72*v0^3*j^3-a0^6+48*x0*j^2*a0
^3-144*x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*x7*j^2*a0^3-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7),
x5 = 1/6*(Y^3+6*x7*j^2)/j^2, v1 = 1/4*(2*Y^2+2*v0*j-a0^2)/j, v5 =
1/2*Y^2/j, a1 = -(-18*Y^2*a0^4+72*Y^2*v0*j*a0^2-72*Y^2*
v0^2*j^2+72*Y*j^2*v0^2*a0+12*Y*a0^5+36*Y*j^2*a0^2*x0-60*Y*j*v0*a0^3-72*Y*j^3*v0*x0-36*Y*j^2*a0^2*x7+72*Y*j^3*
v0*x7+72*x0^2*j^4+72*x7^2*j^4+72*v0^3*j^3-
a0^6+48*x0*j^2*a0^3-144*x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*
x7*j^2*a0^3-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7)*Y/
(72*x0^2*j^4+72*x7^2*j^4+72*v0^3*j^3-a0^6+48*x0*j^2*a0^3-144
*x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*x7*j^2*a0^3-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7),
a5 = Y]
Now, if you really want a formula only involving {j, x0, v0, a0, x7}
you can write out the solutions of the quartic and plug them in,
however I don't recommend this. You might want to refine your
question. Instead of looking for the generic solutions, perhaps you
want values of {j, x0, v0, a0, x7} that produce different types of
solutions ?
Someone who is an expert on solving parametric systems may want to
comment.
pnachtwey
> I took the solution from above and plugged it into the solve command:
>
> Y=RootOf((72*v0*j-36*a0^2)*_Z^4+(-144*x7*j^2-144*v0*a0*j
> +144*x0*j^2+48*a0^3)*_Z^3+(-18*a0^4+72*v0*j*a0^2-72*v0^
> a0^2-48*x0*j^2*a0^3+48*x7*j^2*a0^3+a0^6-6*v0*j*a0^4);
>
> That is, Y is a root of a degree 4 polynomial with coefficients in {j,
> x0, v0, a0, x7}. Then the solutions are:
>
> [x1 =
> *a0+216*x7^2*j^4*a0^3-270*v0^2*j^2*a0^5+216*v0^3*j^5*x7+1728*Y*v0*j^5*x0*x7-1440*Y*v0*j^3*x0*a0^3+1440*Y*v0*j
> ^3*x7*a0^3-1728*Y*v0^2*j^4*a0*x7+1728*Y*v0^2*j^4*x0*a0-864*Y*a0^2*x0*j^4*x7+1080*x0^2*j^4*a0^3+216*Y^3*x0^2*j
> ^4+216*Y^3*x7^2*j^4+1080*Y^3*v0^3*j^3-864*Y*v0^4*j^4-3240*x0^2*j^5*v0*a0-1944*x0^2*j^6*x7+666*Y^3*v0*j*a0^4+
> 864*Y^2*x0*j^4*v0^2-864*Y^2*v0^3*a0*j^3+1152*Y^2*v0^2*a0^3*j^2-216*Y^2*x7*j^2*a0^4-504*Y^2*v0*a0^5*j-864*Y^2*
> x7*j^4*v0^2+216*Y^2*x0*j^2*a0^4+1458*v0*j^3*x7*a0^4-2268*v0^2*j^4*x7*a0^2-111*Y^3*a0^6+72*Y^2*a0^7-6*Y*a0^8+
> 216*x7^3*j^6-5*a0^9+45*v0*j*a0^7+900*v0^3*j^3*a0^3+1620*v0^2*j^4*x0*a0^2-1350*v0*j^3*x0*a0^4+1080*x0^3*j^6-
> \
> 1296*x0*j^4*a0^3*x7+3888*x0*j^5*x7*v0*a0-432*Y^3*x0*j^3*v0*a0+432*Y^3*v0*a0*j^3*x7+864*Y^2*x7*j^3*v0*a0^2-864
> *Y^2*x0*j^3*v0*a0^2+864*Y*v0^3*j^3*a0^2-864*Y*v0*j^5*x0^2+432*Y*a0^2*x7^2*j^4-864*Y*v0*j^5*x7^2+48*Y*v0*j*a0^
> 6+432*Y*a0^2*x0^2*j^4+288*Y*a0^5*x0*j^2-288*Y*a0^5*x7*j^2-288*Y*v0^2*j^2*a0^4-432*Y^3*x0*j^4*x7-1404*Y^3*v0^2
> *j^2*a0^2+144*Y^3*x0*j^2*a0^3-144*Y^3*x7*j^2*a0^3)/j^2/
> (72*x0^2*j^4+72*x7^2*j^4+72*v0^3*j^3-a0^6+48*x0*j^2*a0
> ^3-144*x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*x7*j^2*a0^3-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7),
> x5 = 1/6*(Y^3+6*x7*j^2)/j^2, v1 = 1/4*(2*Y^2+2*v0*j-a0^2)/j, v5 =
> 1/2*Y^2/j, a1 = -(-18*Y^2*a0^4+72*Y^2*v0*j*a0^2-72*Y^2*
> v0*x7+72*x0^2*j^4+72*x7^2*j^4+72*v0^3*j^3-
> a0^6+48*x0*j^2*a0^3-144*x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*
> x7*j^2*a0^3-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7)*Y/
> (72*x0^2*j^4+72*x7^2*j^4+72*v0^3*j^3-a0^6+48*x0*j^2*a0^3-144
> *x0*j^4*x7-36*v0^2*j^2*a0^2+6*v0*j*a0^4-48*x7*j^2*a0^3-144*x0*j^3*v0*a0+144*v0*a0*j^3*x7),
> a5 = Y]
>
> Now, if you really want a formula only involving {j, x0, v0, a0, x7}
> you can write out the solutions of the quartic and plug them in,
> however I don't recommend this. You might want to refine your
> question. Instead of looking for the generic solutions, perhaps you
> want values of {j, x0, v0, a0, x7} that produce different types of
> solutions ?
>
> Someone who is an expert on solving parametric systems may want to
> comment.
from Mathcad, at least the part I can see. The problem I have with
Mathcad is that I can't do anything from there. I get a not responding
error every time I hit a key or try to scroll with the mouse. It does
look like Maple would require some other steps to find the roots just
like Mathcad. Perhaps this is the nature of the answer. I see no
cure for the y=RootOf(poly(z)) problem now. Notice some of the terms
are very simple. It may be possible to solve in terms of them. The
formula for x1 looks nasty.
I tried the same problem with wxMaxima. It didn't even return an
answer.
eq0: x0+v0*(a1-a0)/j+(a0/2)*((a1-a0)/j)^2+(j/6)*((a1-a0)/j)^3=x1;
eq1: v0+a0*(a1-a0)/j+(j/2)*((a1-a0)/j)^2=v1;
eq2: x1+v1*((a1-a5)/j)+(a1/2)*((a1-a5)/j)^2-(j/6)*((a1-a5)/j)^3=x5;
eq3: v1+a1*((a1-a5)/j)+(j/2)*((a1-a5)/j)^2=v5;
eq4: x5+v5*(-a5/j)+(a5/2)*(-a5/j)^2+(j/6)*(-a5/j)^3=x7;
eq5: v5+a5*(-a5/j)+(j/2)*(-a5/j)^2=0;
solve([eq0,eq1,eq2,eq3,eq4,eq5],[x1,v1,a1,x5,v5,a5]);
Thanks, this has been a big help. I will be interesting if anyone
using Mathematica can do better.
Peter Nachtwey
Robert H. Lewis
I used the equations:
x0 + v0*(a1-a0)/j + a0/2*(a1-a0)^2/j^2 + j/6*(a1-a0)^3/j^3 - x1,
v0 + a0*(a1-a0)/j + j/2*(a1-a0)^2/j^2 - v1,
x1 + v1*(a1-a5)/j + a1/2*(a1-a5)^2/j^2 - j/6*(a1-a5)^3/j^3 - x5,
v1 + a1*(a1-a5)/j - j/2*(a1-a5)^2/j^2 - v5,
x5 + v5*(-a5/j) + a5/2*(-a5/j)^2 + j/6*(-a5/j)^3 - x7,
v5 + a5*(-a5/j) + j/2*(-a5/j)^2
where each is set to 0. After simplifying, this is the same as:
a1^3 + 6*v0*j*a1 - 3*a0^2*a1 - 6*j^2*x1 + 6*x0*j^2 - 6*v0*a0*j + 2*a0^3 ,
-2*j*v1 + a1^2 + 2*v0*j - a0^2 ,
a5^3 - 6*j*v1*a5 - 3*a1^2*a5 + 6*j*a1*v1 + 2*a1^3 - 6*j^2*x5 + 6*j^2*x1 ,
- a5^2 + 2*j*v1 + a1^2 - 2*j*v5 ,
a5^3 - 3*j*v5*a5 + 3*j^2*x5 - 3*x7*j^2 ,
- a5^2 + 2*j*v5
There are six variables: x1, v1, a1, x5, v5, a5. The others are parameters. This system is very easy to solve. Solving for a5 yields
72*v0*j*a5^4 - 36*a0^2*a5^4 - 144*x7*j^2*a5^3 + 144*x0*j^2*a5^3 - 144*v0*a0*j*a5^3 + 48*a0^3*a5^3 - 72*v0^2*j^2*a5^2 + 72*v0*a0^2*j*a5^2 - 18*a0^4*a5^2 - 72*x7^2*j^4 + 144*x0*x7*j^4 - 72*x0^2*j^4 - 144*v0*a0*x7*j^3 + 144*x0*v0*a0*j^3 - 72*v0^3*j^3 + 48*a0^3*x7*j^2 - 48*x0*a0^3*j^2 + 36*v0^2*a0^2*j^2 - 6*v0*a0^4*j + a0^6
This takes around 0.01 secs. Solving for x1 is a bit harder. The answer has 244 terms and takes about 0.1 secs. It starts and ends:
483729408*v0^3*j^11*x1^4 - 725594112*v0^2*a0^2*j^10*x1^4 + 362797056*v0*a0^4*j^9*x1^4 ..... 13950*v0*a0^16*j - 775*a0^18
If this is what you had in mind let me know.
This is solved using the method-of-choice for polynomial systems, the Dixon resultant. I used my computer algebra system Fermat. See the web page below for lots of other much harder examples.
Robert H. Lewis
Fordham University
New York
http://home.bway.net/lewis/
pnachtwey
enough. You have to remember that there are few variables here so a
PowerPC with 32 floating point regesters can do calculations with a
minimum of memory accesses. Compilers are very good a combining
common terms. x1 still is nasty. j^11 doesn't look good.
This is what I had in mind though. Now I will look up the Dixon
resultant.
Thanks, I would really like to know the math package you used.
Peter Nachtwey
Robert H. Lewis
> What cas package did you use? The answer for a5
> looks like it is easy
> enough. You have to remember that there are few
> variables here so a
> PowerPC with 32 floating point regesters can do
> calculations with a
> minimum of memory accesses. Compilers are very good
> a combining common terms.
I don't see the point of your comments above.
> x1 still is nasty. j^11 doesn't look
> good.
Here is the complete polynomial. It seems to be irreducible.
483729408*v0^3*j^11*x1^4 - 725594112*v0^2*a0^2*j^10*x1^4 + 362797056*v0*a0^4*j^9*x1^4 - 60466176*a0^6*j^8*x1^4 - 644972544*x7^3*j^12*x1^3
+ 1934917632*x0*x7^2*j^12*x1^3 - 1934917632*x0^2*x7*j^12*x1^3 + 644972544*x0^3*j^12*x1^3 - 1934917632*v0*a0*x7^2*j^11*x1^3 + 3869835264*x0*v0*a0*
x7*j^11*x1^3 - 1451188224*v0^3*x7*j^11*x1^3 - 1934917632*x0^2*v0*a0*j^11*x1^3 - 483729408*x0*v0^3*j^11*x1^3 + 644972544*a0^3*x7^2*j^10*x1^3
- 1289945088*x0*a0^3*x7*j^10*x1^3 + 241864704*v0^2*a0^2*x7*j^10*x1^3 + 644972544*x0^2*a0^3*j^10*x1^3 + 2660511744*x0*v0^2*a0^2*j^10*x1^3
+ 483729408*v0^4*a0*j^10*x1^3 + 201553920*v0*a0^4*x7*j^9*x1^3 - 1652742144*x0*v0*a0^4*j^9*x1^3 - 1531809792*v0^3*a0^3*j^9*x1^3 - 33592320*a0^6*
x7*j^8*x1^3 + 275457024*x0*a0^6*j^8*x1^3 + 1249634304*v0^2*a0^5*j^8*x1^3 - 396389376*v0*a0^7*j^7*x1^3 + 44043264*a0^9*j^6*x1^3 + 161243136*x7^4*
j^12*x1^2 + 1289945088*x0*x7^3*j^12*x1^2 - 4837294080*x0^2*x7^2*j^12*x1^2 + 5159780352*x0^3*x7*j^12*x1^2 - 1773674496*x0^4*j^12*x1^2 - 1289945088*v0*
a0*x7^3*j^11*x1^2 + 9674588160*x0*v0*a0*x7^2*j^11*x1^2 + 120932352*v0^3*x7^2*j^11*x1^2 - 15479341056*x0^2*v0*a0*x7*j^11*x1^2 + 4111699968*x0*v0^3*
x7*j^11*x1^2 + 7094697984*x0^3*v0*a0*j^11*x1^2 - 1330255872*x0^2*v0^3*j^11*x1^2 + 429981696*a0^3*x7^3*j^10*x1^2 - 3224862720*x0*a0^3*x7^2*j^10*x1^2
- 5018692608*v0^2*a0^2*x7^2*j^10*x1^2 + 5159780352*x0^2*a0^3*x7*j^10*x1^2 + 9311791104*x0*v0^2*a0^2*x7*j^10*x1^2 - 4111699968*v0^4*a0*x7*j^10*x1^2
- 2364899328*x0^3*a0^3*j^10*x1^2 - 8646663168*x0^2*v0^2*a0^2*j^10*x1^2 + 2660511744*x0*v0^4*a0*j^10*x1^2 + 53747712*v0^6*j^10*x1^2 + 3315561984*v0*
a0^4*x7^2*j^9*x1^2 - 7235785728*x0*v0*a0^4*x7*j^9*x1^2 + 2378336256*v0^3*a0^3*x7*j^9*x1^2 + 6097006080*x0^2*v0*a0^4*j^9*x1^2 + 2217093120*x0*
v0^3*a0^3*j^9*x1^2 - 1491499008*v0^5*a0^2*j^9*x1^2 - 552593664*a0^6*x7^2*j^8*x1^2 + 1205964288*x0*a0^6*x7*j^8*x1^2 + 20155392*v0^2*a0^5*x7*j^8*x1^2
- 1016167680*x0^2*a0^6*j^8*x1^2 - 3769058304*x0*v0^2*a0^5*j^8*x1^2 + 1310100480*v0^4*a0^4*j^8*x1^2 - 178039296*v0*a0^7*x7*j^7*x1^2 + 1367207424*
x0*v0*a0^7*j^7*x1^2 - 245223936*v0^3*a0^6*j^7*x1^2 + 19782144*a0^9*x7*j^6*x1^2 - 151911936*x0*a0^9*j^6*x1^2 - 78941952*v0^2*a0^8*j^6*x1^2
+ 30979584*v0*a0^10*j^5*x1^2 - 2581632*a0^12*j^4*x1^2 - 13436928*x7^5*j^12*x1 - 255301632*x0*x7^4*j^12*x1 - 779341824*x0^2*x7^3*j^12*x1 + 4004204544*
x0^3*x7^2*j^12*x1 - 4581992448*x0^4*x7*j^12*x1 + 1625868288*x0^5*j^12*x1 + 255301632*v0*a0*x7^4*j^11*x1 + 1558683648*x0*v0*a0*x7^3*j^11*x1
- 600182784*v0^3*x7^3*j^11*x1 - 12012613632*x0^2*v0*a0*x7^2*j^11*x1 + 1558683648*x0*v0^3*x7^2*j^11*x1 + 18327969792*x0^3*v0*a0*x7*j^11*x1
- 5670383616*x0^2*v0^3*x7*j^11*x1 - 8129341440*x0^4*v0*a0*j^11*x1 + 2776965120*x0^3*v0^3*j^11*x1 - 85100544*a0^3*x7^4*j^10*x1 - 519561216*x0*a0^3*
x7^3*j^10*x1 + 120932352*v0^2*a0^2*x7^3*j^10*x1 + 4004204544*x0^2*a0^3*x7^2*j^10*x1 + 9674588160*x0*v0^2*a0^2*x7^2*j^10*x1 - 1558683648*v0^4*a0*x7^2*
j^10*x1 - 6109323264*x0^3*a0^3*x7*j^10*x1 - 18986379264*x0^2*v0^2*a0^2*x7*j^10*x1 + 11340767232*x0*v0^4*a0*x7*j^10*x1 - 1303382016*v0^6*x7*j^10*
x1 + 2709780480*x0^4*a0^3*j^10*x1 + 12093235200*x0^3*v0^2*a0^2*j^10*x1 - 8330895360*x0^2*v0^4*a0*j^10*x1 + 1195886592*x0*v0^6*j^10*x1
+ 69424128*v0*a0^4*x7^3*j^9*x1 - 6839396352*x0*v0*a0^4*x7^2*j^9*x1 - 1146617856*v0^3*a0^3*x7^2*j^9*x1 + 14075182080*x0^2*v0*a0^4*x7*j^9*x1
- 2463436800*x0*v0^3*a0^3*x7*j^9*x1 - 1760237568*v0^5*a0^2*x7*j^9*x1 - 8756398080*x0^3*v0*a0^4*j^9*x1 - 985374720*x0^2*v0^3*a0^3*j^9*x1
+ 4743235584*x0*v0^5*a0^2*j^9*x1 - 1195886592*v0^7*a0*j^9*x1 - 11570688*a0^6*x7^3*j^8*x1 + 1139899392*x0*a0^6*x7^2*j^8*x1 + 2055849984*v0^2*a0^5*x7^2*j^8*x1 - 2345863680*x0^2*a0^6*x7*j^8*x1 - 4152010752*x0*v0^2*a0^5*x7*j^8*x1 + 2816156160*v0^4*a0^4*x7*j^8*x1 + 1459399680*x0^3*a0^6*
j^8*x1 + 5845063680*x0^2*v0^2*a0^5*j^8*x1 - 5436357120*x0*v0^4*a0^4*j^8*x1 + 1209323520*v0^6*a0^3*j^8*x1 - 750228480*v0*a0^7*x7^2*j^7*x1
+ 1856535552*x0*v0*a0^7*x7*j^7*x1 - 1185435648*v0^3*a0^6*x7*j^7*x1 - 2295475200*x0^2*v0*a0^7*j^7*x1 + 1675883520*x0*v0^3*a0^6*j^7*x1
- 363916800*v0^5*a0^5*j^7*x1 + 83358720*a0^9*x7^2*j^6*x1 - 206281728*x0*a0^9*x7*j^6*x1 + 212471424*v0^2*a0^8*x7*j^6*x1 + 255052800*x0^2*a0^9*j^6*x1 - 54587520*x0*v0^2*a0^8*j^6*x1 + 20528640*v0^4*a0^7*j^6*x1 - 21866112*v0*a0^10*x7*j^5*x1 - 40093056*x0*v0*a0^10*j^5*x1 - 3058560*v0^3*a0^9*
j^5*x1 + 1822176*a0^12*x7*j^4*x1 + 3341088*x0*a0^12*j^4*x1 + 4478976*v0^2*a0^11*j^4*x1 - 946080*v0*a0^13*j^3*x1 + 63072*a0^15*j^2*x1 + 373248*x7^6*j^12 + 11197440*x0*x7^5*j^12 + 99657216*x0^2*x7^4*j^12 + 126904320*x0^3*x7^3*j^12 - 1096229376*x0^4*x7^2*j^12 + 1354890240*x0^5*x7*
j^12 - 496793088*x0^6*j^12 - 11197440*v0*a0*x7^5*j^11 - 199314432*x0*v0*a0*x7^4*j^11 + 164602368*v0^3*x7^4*j^11 - 380712960*x0^2*v0*a0*x7^3*j^11
- 58226688*x0*v0^3*x7^3*j^11 + 4384917504*x0^3*v0*a0*x7^2*j^11 - 692001792*x0^2*v0^3*x7^2*j^11 - 6774451200*x0^4*v0*a0*x7*j^11 + 2351462400*x0^3*v0^3*x7*j^11 + 2980758528*x0^5*v0*a0*j^11 - 1282106880*x0^4*v0^3*j^11 + 3732480*a0^3*x7^5*j^10 + 66438144*x0*a0^3*x7^4*j^10 - 147246336*v0^2*a0^2*
x7^4*j^10 + 126904320*x0^2*a0^3*x7^3*j^10 + 468052992*x0*v0^2*a0^2*x7^3*j^10 + 58226688*v0^4*a0*x7^3*j^10 - 1461639168*x0^3*a0^3*x7^2*j^10
- 5539373568*x0^2*v0^2*a0^2*x7^2*j^10 + 1384003584*x0*v0^4*a0*x7^2*j^10 + 148925952*v0^6*x7^2*j^10 + 2258150400*x0^4*a0^3*x7*j^10 + 10021708800*
x0^3*v0^2*a0^2*x7*j^10 - 7054387200*x0^2*v0^4*a0*x7*j^10 + 1005530112*x0*v0^6*x7*j^10 - 993586176*x0^5*a0^3*j^10 - 5528736000*x0^4*v0^2*a0^2*j^10
+ 5128427520*x0^3*v0^4*a0*j^10 - 1100708352*x0^2*v0^6*j^10 + 57013632*v0*a0^4*x7^4*j^9 - 297478656*x0*v0*a0^4*x7^3*j^9 - 233653248*v0^3*a0^3*x7^3* j^9 + 3865916160*x0^2*v0*a0^4*x7^2*j^9 + 1847577600*x0*v0^3*a0^3*x7^2*j^9 - 1138779648*v0^5*a0^2*x7^2*j^9 - 7269004800*x0^3*v0*a0^4*x7*j^9
- 615859200*x0^2*v0^3*a0^3*x7*j^9 + 4037796864*x0*v0^5*a0^2*x7*j^9 - 1005530112*v0^7*a0*x7*j^9 + 4006350720*x0^4*v0*a0^4*j^9 + 533744640*x0^3*
v0^3*a0^3*j^9 - 4390516224*x0^2*v0^5*a0^2*j^9 + 2201416704*x0*v0^7*a0*j^9 - 313901568*v0^9*j^9 - 9502272*a0^6*x7^4*j^8 + 49579776*x0*a0^6*x7^3*
j^8 + 199687680*v0^2*a0^5*x7^3*j^8 - 644319360*x0^2*a0^6*x7^2*j^8 - 2654913024*x0*v0^2*a0^5*x7^2*j^8 + 961580160*v0^4*a0^4*x7^2*j^8 + 1211500800*
x0^3*a0^6*x7*j^8 + 4730918400*x0^2*v0^2*a0^5*x7*j^8 - 4739316480*x0*v0^4*a0^4*x7*j^8 + 1000304640*v0^6*a0^3*x7*j^8 - 667725120*x0^4*a0^6*j^8 - 3525327360*x0^3*v0^2*a0^5*j^8 + 5087836800*x0^2*v0^4*a0^4*j^8 - 2209628160*x0*v0^6*a0^3*j^8 + 311848704*v0^8*a0^2*j^8 - 64136448*v0*a0^7*x7^3*
j^7 + 942637824*x0*v0*a0^7*x7^2*j^7 -198567936*v0^3*a0^6*x7^2*j^7 - 1870905600*x0^2*v0*a0^7*x7*j^7 + 1582571520*x0*v0^3*a0^6*x7*j^7 - 252502272*v0^5*a0^5*x7*j^7 + 1388793600*x0^3*v0*a0^7*j^7 - 1629227520*x0^2*v0^3*a0^6*j^7 + 616419072*x0*v0^5*a0^5*j^7 - 71290368*v0^7*a0^4*j^7
+ 7126272*a0^9*x7^3*j^6 - 104737536*x0*a0^9*x7^2*j^6 - 43366752*v0^2*a0^8*x7^2*j^6 + 207878400*x0^2*a0^9*x7*j^6 - 125737920*x0*v0^2*a0^8*x7*j^6
- 45722880*v0^4*a0^7*x7*j^6 - 154310400*x0^3*a0^9*j^6 + 90162720*x0^2*v0^2*a0^8*j^6 + 25194240*x0*v0^4*a0^7*j^6 - 19564416*v0^6*a0^6*j^6
+ 19147104*v0*a0^10*x7^2*j^5 - 16428096*x0*v0*a0^10*x7*j^5 + 34292160*v0^3*a0^9*x7*j^5 + 28260576*x0^2*v0*a0^10*j^5 - 31233600*x0*v0^3*a0^9*j^5
+ 11524032*v0^5*a0^8*j^5 - 1595592*a0^12*x7^2*j^4 + 1369008*x0*a0^12*x7*j^4 - 7858944*v0^2*a0^11*x7*j^4 - 2355048*x0^2*a0^12*j^4 + 3379968*x0*v0^2*
a0^11*j^4 - 2638656*v0^4*a0^10*j^4 + 1103760*v0*a0^13*x7*j^3 - 157680*x0*v0*a0^13*j^3 + 597888*v0^3*a0^12*j^3 - 73584*a0^15*x7*j^2 + 10512*x0*a0^15*
j^2 - 116856*v0^2*a0^14*j^2 + 13950*v0*a0^16*j - 775*a0^18
> This is what I had in mind though. Now I will look
> up the Dixon resultant.
>
> Thanks, I would really like to know the math package
> you used.
>
> Peter Nachtwey
>
I used Fermat. I wrote it. You can find it at the link below. But I don't that's the real point here.
Robert H. Lewis
Fordham University
Peter Nachtwey
"Robert H. Lewis" <rle...@fordham.edu> wrote in message
news:6773056.12071059050...@nitrogen.mathforum.org...
Ouch, that isn't going to happen real time.
> I used Fermat. I wrote it. You can find it at the link below. But I don't
> that's the real point here.
I suspected that you used something other than Maple or Mathematica.. I
usually search the internet for the names of people I respond to or respond
to me. I found that you have a .pdf explaining this very topic so I am not
surprised that you have written your own symbolic solver.
So far I find this depressing. I think I can solve this iteratively faster
than doing all those calculations. I know that a1 must be between 0 and the
requested acceleration provide by the command to the motion controller. The
problem is that a motion controller must be able to solve these problems in
a few microseconds. Clearly this will not happen. :( What gets me is that
the solution is so simple if the initial velocity and acceleration are 0.
I will look at Fermat. I noticed you have a Mac version.
Thanks again.
Peter Nachtwey
Bob Walton
> I want to generate a simple 3rd order motion profile that will go from
> an initial PVA ( position, velocity and acceleration ) to a final PVA
> where the final V=0 and A=0. Sound easy, I have 6 equations and 6
> unknows to solve for. There are three motion segments. The first is
> from point 0 to point 1 where the jerk is postive and the target is
> accelerating. The second is from point 1 to point 5. ( i skipped
> some points ) where the jerk is negative and a5 is at peak
> deceleration. The final segment goes from point 5 to point 7. Duing
> this segment the jerk is postive again and the deceleration is
> decreasing until the target stops at x7.
>
> Here is a link to my Mathcad work sheet.
> ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Seg137.pdf
> It looks like it should be simple. Mathcad does come up with an answer
> but after that it doesn't respond to the keyboard or mouse very well
> or at all. I keep getting a Not Responding errorr. I can't cut and
> paste the resulting RootOf equation so I can't solve for the roots and
> do the manual operations. BTW, I hate having to cut and paste
> equations or intermediate results when I change a parameter.
>
> Will Mathematica, Maple or other CAS solve this, symbolically,
> easily, quickly? Can it produce equations in text form so I can cut
> and paste it into my C program? At this time I must retype all the
> Mathcad formulas into C. This is tedious and error prone.
>
Well, I used Maxima as a scratchpad for a manual solution. The four
variables v1, v5, x1 and x5 are easily eliminated automatically. Of the
two resulting equations, the second is linear in a1, and easily solved.
The result is substituted into the first equation, resulting in a
quartic in a5. This is best solved by a good numerical solver, such as
Matlab. Here is the Maxima and Matlab code I used:
/*Maxima code*/
eq0:x0+v0*(a1-a0)/j+(a0/2)*((a1-a0)/j)^2+(j/6)*
((a1-a0)/j)^3-x1;
eq1:v0+a0*(a1-a0)/j+(j/2)*((a1-a0)/j)^2-v1;
eq2:x1+v1*((a1-a5)/j)+(a1/2)*((a1-a5)/j)^2-
(j/6)*((a1-a5)/j)^3-x5;
eq3:v1+a1*((a1-a5)/j)-(j/2)*((a1-a5)/j)^2-v5;
eq4:x5+v5*(-a5/j)+(a5/2)*(-a5/j)^2+(j/6)*
(-a5/j)^3-x7;
eq5:v5+a5*(-a5/j)+(j/2)*(-a5/j)^2;
q6:eliminate([eq0,eq1,eq2,eq3,eq4,eq5],
[v5,v1,x5,x1]);
q7:solve(q6[2],a1);
q8:subst(q7,q6[1]);
q9:num(ratsimp(q8)); /*we only care about the numerator*/
c4:ratcoef(q9,a5,4); /*coefficients of quartic*/
c3:ratcoef(q9,a5,3);
c2:ratcoef(q9,a5,2);
c1:ratcoef(q9,a5,1);
c0:ratcoef(q9,a5,0);
%Matlab code
x0=10; %garbage set of input data for testing
v0=10;
j=10;
a0=10;
x7=100;
c4=288*j^3*v0-144*a0^2*j^2; %4th order term of quadratic
c3=-576*j^4*x7+576*j^4*x0-576*a0*j^3*v0+192*a0^3*j^2; %3rd order
c2=-288*j^4*v0^2+288*a0^2*j^3*v0-72*a0^4*j^2; %2nd order
c1=0; %1st order
c0=-288*j^6*x7^2+(576*j^6*x0-576*a0*j^5*v0+192*a0^3*j^4)...
*x7-288*j^6*x0^2+(576*a0*j^5*v0-192*a0^3*j^4)*x0-...
288*j^5*v0^3+144*a0^2*j^4*v0^2-24*a0^4*j^3*v0+4*a0^6*j^2;
a5=roots([c4 c3 c2 c1 c0]) %solve for four roots,
%carry all through
a1=(6*j^2*x7-6*j^2*x0+6*a0*j*v0+6*a5.^3-2*a0^3)./...
(6*j*v0+6*a5.^2-3*a0^2) %q7 from Maxima results
v5=a5.^2/(2*j) %solve(eq5,v5)
x5=(3*j^2*x7+3*a5*j.*v5-a5.^3)/(3*j^2) %solve(eq4,x5)
x1=(6*j^2*x0+(6*a1-6*a0)*j*v0+a1.^3-3*a0^2*a1+2*a0^3)...
/(6*j^2) %solve(eq0,x1)
v1=(2*j*v0+a1.^2-a0^2)/(2*j) %solve(eq1,v1)
%demonstrate that the set of equations is in
%fact solved by all roots
eq0=-(6*j^2*x1-6*j^2*x0+(6*a0-6*a1)*j*v0-a1.^3+...
3*a0^2*a1-2*a0^3)/(6*j^2)
eq1=-(2*j*v1-2*j*v0-a1.^2+a0^2)/(2*j)
eq2=-(6*j^2*x5-6*j^2*x1+(6*a5-6*a1)*j.*v1-a5.^3+...
3*a1.^2.*a5-2*a1.^3)/(6*j^2)
eq3=-(2*j*v5-2*j*v1+a5.^2-a1.^2)/(2*j)
eq4=-(3*j^2*x7-3*j^2*x5+3*a5*j.*v5-a5.^3)/(3*j^2)
eq5=(2*j*v5-a5.^2)/(2*j)
In this case, only the first root was real and positive, so
for this case it is clear which root is desired.
One may use Maxima's "grind" function to generate code that can
be easily copy/pasted into an editor for use with other software,
such as Matlab in my example above.
...
> Peter Nachtwey
--
Bob Walton
Email: http://bwalton.com/cgi-bin/emailbob.pl
pnachtwey
> The four
> variables v1, v5, x1 and x5 are easily eliminated automatically.
Yes, I could have done that manually in my example using substitution.
> Of the
> two resulting equations, the second is linear in a1, and easily solved.
> The result is substituted into the first equation, resulting in a
> quartic in a5. This is best solved by a good numerical solver, such as
> Matlab. Here is the Maxima and Matlab code I used:
>
> /*Maxima code*/
> eq0:x0+v0*(a1-a0)/j+(a0/2)*((a1-a0)/j)^2+(j/6)*
> ((a1-a0)/j)^3-x1;
> eq1:v0+a0*(a1-a0)/j+(j/2)*((a1-a0)/j)^2-v1;
> eq2:x1+v1*((a1-a5)/j)+(a1/2)*((a1-a5)/j)^2-
> (j/6)*((a1-a5)/j)^3-x5;
> eq3:v1+a1*((a1-a5)/j)-(j/2)*((a1-a5)/j)^2-v5;
> eq4:x5+v5*(-a5/j)+(a5/2)*(-a5/j)^2+(j/6)*
> (-a5/j)^3-x7;
> eq5:v5+a5*(-a5/j)+(j/2)*(-a5/j)^2;
> q6:eliminate([eq0,eq1,eq2,eq3,eq4,eq5],
> [v5,v1,x5,x1]);
> q7:solve(q6[2],a1);
I wonder if this trick will help with the Mathcad
> q8:subst(q7,q6[1]);
and this too
> q9:num(ratsimp(q8)); /*we only care about the numerator*/
> c4:ratcoef(q9,a5,4); /*coefficients of quartic*/
> c3:ratcoef(q9,a5,3);
> c2:ratcoef(q9,a5,2);
> c1:ratcoef(q9,a5,1);
> c0:ratcoef(q9,a5,0);
This part is obvious and you are right about the real root. I used
these equations to generate the coefficients for the polyroots
function in Mathcad. The real root is the one I want. I then
substituted back into a1, and the other variables and plotted the
results. It worked!!!
I don't have Matlab buf I do have Scilab. It looks like between
Maxima and Scilab one can do a lot.
>
> One may use Maxima's "grind" function to generate code that can
> be easily copy/pasted into an editor for use with other software,
> such as Matlab in my example above.
This is good. I liked the technique. I wonder if Mathcad can solve
this if I manually do the elimination and just substite for a1. I
will give it a try. Meanwhile I must also find an efficient way of
solving the roots or I must try a different motion profile all
together.
Thanks,
Peter Nachtwey
pnachtwey
problems right away. Mathcad does not have an eliminate function so I
tried just using substitution. I simplified the 6 equations into just
two One is for the positions and the other is for the velocities as a
function of x0,v0,a0,a1,a5,x7. I verified that my formulas for
positions and velocity were still correct using the numbers I
generated with your formulas. However, my formula for position and
velocity are NOT the same as the two formulas your example provided
for q6. The equation for the velocity is:
(2*v0*j-a0^2+2*a1^2-2*a5^2)/(2*j)
and the position formula is:
(3*(x0-x7)*j^2+(-6*v0*a5-3*v0*a0+6*v0*a1)*j
+3*a0^2*a5+a0^3+3*a1^3-6*a1^2*a5-3*a.0^2*a1+3*a5^3)/(3*j^2)
These equations are equal to 0 (close, because I am cutting and
pasting numbers) when using the accelerations generated by your
formulas. You can see there are some similarities between my velocity
formula and q6[1] but there are not the same. q6[1] doesn't even
generate a value of 0 when real numbers are used for v0,a1,a5 and j.
There must be some trick to the eliminate function where is does
something other than just making substitutions to remove variables.
Also, when I solve for a1 and a5 I get a mess which is much more
complicated than the coefficients your method generated.
I read the documentation in the Maxima help. It sounds like it is
just doing simple substitutions but that can't be.
Does any body have a explanation for what the eliminate function
really does?
Peter Nachtwey
Bob Walton
> I tried to copy your technique and apply it to Mathcad. I ran into
> problems right away. Mathcad does not have an eliminate function so I
> tried just using substitution. I simplified the 6 equations into just
Maxima is free -- sourceforge.net
> two One is for the positions and the other is for the velocities as a
> function of x0,v0,a0,a1,a5,x7. I verified that my formulas for
> positions and velocity were still correct using the numbers I
> generated with your formulas. However, my formula for position and
> velocity are NOT the same as the two formulas your example provided
> for q6. The equation for the velocity is:
>
> (2*v0*j-a0^2+2*a1^2-2*a5^2)/(2*j)
>
> and the position formula is:
>
> (3*(x0-x7)*j^2+(-6*v0*a5-3*v0*a0+6*v0*a1)*j
> +3*a0^2*a5+a0^3+3*a1^3-6*a1^2*a5-3*a.0^2*a1+3*a5^3)/(3*j^2)
>
> These equations are equal to 0 (close, because I am cutting and
> pasting numbers) when using the accelerations generated by your
> formulas. You can see there are some similarities between my velocity
> formula and q6[1] but there are not the same. q6[1] doesn't even
> generate a value of 0 when real numbers are used for v0,a1,a5 and j.
> There must be some trick to the eliminate function where is does
> something other than just making substitutions to remove variables.
>
> Also, when I solve for a1 and a5 I get a mess which is much more
> complicated than the coefficients your method generated.
Yes, I got an equation (my q6[2]) which was linear in a1, which means
that a1 was easily solved for and substituted into q6[1], resulting in a
quartic in a5 alone -- an equation like that is the goal.
q6[1] is
-4*j^2*(2*j*v0-2*a5^2+2*a1^2-a0^2),
q6[2] is
-576*j^9*(6*j^2*x7-6*j^2*x0+j*(6*a0*v0-6*a1*v0)+
6*a5^3-6*a1*a5^2+3*a0^2*a1-2*a0^3)
Solving q6[2] for a1 gives:
a1 = (6*j^2*x7-6*j^2*x0+6*a0*j*v0+6*a5^3-2*a0^3)/
(6*j*v0+6*a5^2-3*a0^2)
and the numerator of the result of substituting that into q6[1] is:
-288*j^6*x7^2-(-576*j^6*x0+576*a0*j^5*v0+
(576*a5^3-192*a0^3)*j^4)*x7-288*j^6*x0^2-
((192*a0^3-576*a5^3)*j^4-576*a0*j^5*v0)*x0
-288*j^5*v0^3-(288*a5^2-144*a0^2)*j^4*v0^2-
(-288*a5^4+576*a0*a5^3-288*a0^2*a5^2+24*a0^4)*j^3*v0
-(144*a0^2*a5^4-192*a0^3*a5^3+72*a0^4*a5^2-4*a0^6)*j^2
which is a quartic in a5 with only inputs appearing in the coefficients.
>
> I read the documentation in the Maxima help. It sounds like it is
> just doing simple substitutions but that can't be.
>
> Does any body have a explanation for what the eliminate function
> really does?
>
Well, there is a very large set of possible linear combinations of
equations that could result from various attempts at elimination. For
example, if Maxima's eliminate() function is asked to also do a1 after
the other four, it generates an equation that has j^20 through j^24 in
it, and also some quite large coefficients. One can divide by j^20 and
one of the coefficients to get something more reasonable. It's not
wrong, it's just a pain -- although some of the more bizarre sets of
equations may give numerical problems during evaluation. I note that
the final equation for a5 between the two cases (my original and the one
with a1 added last to eliminate()'s variable list) simply differ by a
factor of 331776*j^18.
The order of elimination matters also -- asking for the more complicated
ones first will cause trouble -- if a1 is asked for first, eliminate()
eventually bombs out on lack of memory (on a 2 Gb machine!).
I know nothing about the internals of eliminate(), but Maxima is
open-source, so one could find out exactly what it does.
I have found it generally best to avoid having a CAS solve cubic or
quartic polynomials (or even quadratics) symbolically -- it will do it,
but the results are not useful (usually they are numerically
ill-conditioned besides being super humongo) -- not nearly as useful as
taking the symbolic coefficients of the polynomial from the CAS and
having a good numerical package find the roots. There are zillions of
good root-finders, and any of them could be used for that task. I
suppose Maxima might even have a numeric root-finder that could be used
(find_root, maybe, or realroots or allroots?), but I have learned over
the years to use tools for what they are best at. Matlab has a
world-class univariate polynomial root-finder.
pnachtwey
function.
>
> Yes, I got an equation (my q6[2]) which was linear in a1, which means
> that a1 was easily solved for and substituted into q6[1], resulting in a
> quartic in a5 alone -- an equation like that is the goal.
>
> q6[1] is
>
> -4*j^2*(2*j*v0-2*a5^2+2*a1^2-a0^2),
>
> q6[2] is
>
> -576*j^9*(6*j^2*x7-6*j^2*x0+j*(6*a0*v0-6*a1*v0)+
> 6*a5^3-6*a1*a5^2+3*a0^2*a1-2*a0^3)
>
> Solving q6[2] for a1 gives:
> a1 = (6*j^2*x7-6*j^2*x0+6*a0*j*v0+6*a5^3-2*a0^3)/
> (6*j*v0+6*a5^2-3*a0^2)
>
> and the numerator of the result of substituting that into q6[1] is:
>
> -288*j^6*x7^2-(-576*j^6*x0+576*a0*j^5*v0+
> (576*a5^3-192*a0^3)*j^4)*x7-288*j^6*x0^2-
> ((192*a0^3-576*a5^3)*j^4-576*a0*j^5*v0)*x0
> -288*j^5*v0^3-(288*a5^2-144*a0^2)*j^4*v0^2-
> (-288*a5^4+576*a0*a5^3-288*a0^2*a5^2+24*a0^4)*j^3*v0
> -(144*a0^2*a5^4-192*a0^3*a5^3+72*a0^4*a5^2-4*a0^6)*j^2
>
> which is a quartic in a5 with only inputs appearing in the coefficients.
>
eliminate x1,v1,x5,v5? There is more to that just substitution. See
my link below where I have your solution plotted out and my attempts
to use Mathcad's substitution to do the elimination.
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Seg137%20Bob%20Walton%20and%20Maxima.pdf
> > Does any body have a explanation for what the eliminate function
> > really does?
>
> Well, there is a very large set of possible linear combinations of
> equations that could result from various attempts at elimination. For
> example, if Maxima's eliminate() function is asked to also do a1 after
> the other four, it generates an equation that has j^20 through j^24 in
> it, and also some quite large coefficients. One can divide by j^20 and
> one of the coefficients to get something more reasonable. It's not
> wrong, it's just a pain -- although some of the more bizarre sets of
> equations may give numerical problems during evaluation. I note that
> the final equation for a5 between the two cases (my original and the one
> with a1 added last to eliminate()'s variable list) simply differ by a
> factor of 331776*j^18.
I will have to think about this. What I did using Mathcad's
substitution
sense to me. I don't see how any other answer can be right because
the
resulting equations are eqaul to 0. Yours are not. Yet using your
equations results in finding the right roots.
>
> The order of elimination matters also -- asking for the more complicated
> ones first will cause trouble -- if a1 is asked for first, eliminate()
> eventually bombs out on lack of memory (on a 2 Gb machine!).
I can see that IF if had a choice between a linear and none linear
equation but in my case none of the equations are linear in a1 or a5.
>
> I know nothing about the internals of eliminate(), but Maxima is
> open-source, so one could find out exactly what it does.
>
I was just hoping you are someone would know how it got solutions that
are different from mine.
> I have found it generally best to avoid having a CAS solve cubic or
> quartic polynomials (or even quadratics) symbolically -- it will do it,
> but the results are not useful (usually they are numerically
> ill-conditioned besides being super humongo) -- not nearly as useful as
> taking the symbolic coefficients of the polynomial from the CAS and
> having a good numerical package find the roots.
Actually, I was going to test Robert H Lewis's solution next.
His solution is the only one that doesn't require finding roots.
The problem is that I am not familiar with Fermat and Rober H Lewis
didn't
post the script like you did. Also, I liked the idea of solving for
the
acceleratons first and then substituting. If Fermat can calculate any
of the accelerations without finding roots then I may have a answer.
I have other ways of solving this problem that have either more
constraints
or have more segments. This problem is one that just looks easy at
first
and turns out to be a difficult problem. I have similar simple
problems to solve so the answer to this specific problem is NOT as
important as finding a general solution to this class of problems.
> There are zillions of
> good root-finders, and any of them could be used for that task. I
> suppose Maxima might even have a numeric root-finder that could be used
> (find_root, maybe, or realroots or allroots?), but I have learned over
> the years to use tools for what they are best at. Matlab has a
> world-class univariate polynomial root-finder.
The problem with finding roots is that it must be done realtime in an
embedded product. The good part is that it is easy to estimate where
the root should be so only a few iterations should be required but I
like the certainty of solving an equation better.
Thanks again
Peter Nachtwey
Bob Walton
OK, here is a solution involving just substitution. The numeric answers
are the same, but the coefficients of the quartic are radically
different -- it even solves for a different variable (a1).
Maxima code:
eq0:x0+v0*(a1-a0)/j+(a0/2)*((a1-a0)/j)^2+(j/6)*
((a1-a0)/j)^3-x1;
eq1:v0+a0*(a1-a0)/j+(j/2)*((a1-a0)/j)^2-v1;
eq2:x1+v1*((a1-a5)/j)+(a1/2)*((a1-a5)/j)^2-(j/6)*
((a1-a5)/j)^3-x5;
eq3:v1+a1*((a1-a5)/j)-(j/2)*((a1-a5)/j)^2-v5;
eq4:x5+v5*(-a5/j)+(a5/2)*(-a5/j)^2+(j/6)*(-a5/j)^3-x7;
eq5:v5+a5*(-a5/j)+(j/2)*(-a5/j)^2;
q1:solve(eq5,v5);
q2:solve(eq4,x5); /*contains v5*/
q3:solve(eq0,x1);
q4:solve(eq1,v1);
q5:subst(q1,q2); /*get rid of v5: x5=f(a5)*/
q6:subst(q1,eq3);
q6:subst(q4,q6); /*a1^2 and a5^2*/
q7:subst(q3,eq2); /*remove x1*/
q7:subst(q5,q7); /*remove x5*/
q7:subst(q4,q7); /*remove v1*/
q8:solve(q7,a5); /*it was linear in a5, so solve for a5*/
q9:num(ratsimp(subst(q8,q6))); /*a quartic in a1*/
c4:grind(ratsimp(ratcoef(q9,a1,4)));
c3:grind(ratsimp(ratcoef(q9,a1,3)));
c2:grind(ratsimp(ratcoef(q9,a1,2)));
c1:grind(ratsimp(ratcoef(q9,a1,1)));
c0:grind(ratsimp(ratcoef(q9,a1,0)));
grind(ratsimp(q8));
grind(ratsimp(q5));
grind(ratsimp(q1));
grind(ratsimp(q3));
grind(ratsimp(q4));
grind(ratsimp(eq0));
grind(ratsimp(eq1));
grind(ratsimp(eq2));
grind(ratsimp(eq3));
grind(ratsimp(eq4));
grind(ratsimp(eq5));
Matlab code:
x0=10;
v0=10;
j=10;
a0=10;
x7=100;
c4=36*a0^2-72*j*v0;
c3=144*j^2*x7-144*j^2*x0+144*a0*j*v0-48*a0^3;
c2=-72*j^2*v0^2+72*a0^2*j*v0-18*a0^4;
c1=(288*j^3*v0-144*a0^2*j^2)*x7+(144*a0^2*j^2-288*j^3*v0)*...
x0+288*a0*j^2*v0^2-240*a0^3*j*v0+48*a0^5;
c0=-72*j^4*x7^2+(144*j^4*x0-144*a0*j^3*v0+48*a0^3*j^2)*...
x7-72*j^4*x0^2+(144*a0*j^3*v0-48*a0^3*j^2)*x0+72*j^3*v0^3-...
180*a0^2*j^2*v0^2+102*a0^4*j*v0-17*a0^6;
a1=roots([c4 c3 c2 c1 c0])
a5=-(6*j^2*x7-6*j^2*x0+(6*a0-12*a1)*j*v0-6*a1.^3+6*a0^2*a1-...
2*a0^3)./(6*j*v0+6*a1.^2-3*a0^2)
x5=(6*j^2*x7+a5.^3)/(6*j^2)
v5=a5.^2/(2*j)
x1=(6*j^2*x0+(6*a1-6*a0)*j*v0+a1.^3-3*a0^2*a1+2*a0^3)/(6*j^2)
v1=(2*j*v0+a1.^2-a0^2)/(2*j)
%prove the solutions
eq0=-(6*j^2*x1-6*j^2*x0+(6*a0-6*a1)*j*v0-a1.^3+3*a0^2*a1-...
2*a0^3)/(6*j^2)
eq1=-(2*j*v1-2*j*v0-a1.^2+a0^2)/(2*j)
eq2=-(6*j^2*x5-6*j^2*x1+(6*a5-6*a1)*j.*v1-a5.^3+3*a1.^2.*a5-...
2*a1.^3)/(6*j^2)
eq3=-(2*j*v5-2*j*v1+a5.^2-a1.^2)/(2*j)
eq4=-(3*j^2*x7-3*j^2*x5+3*a5*j.*v5-a5.^3)/(3*j^2)
eq5=(2*j*v5-a5.^2)/(2*j)
>
>> The order of elimination matters also -- asking for the more complicated
>> ones first will cause trouble -- if a1 is asked for first, eliminate()
>> eventually bombs out on lack of memory (on a 2 Gb machine!).
>
> I can see that IF if had a choice between a linear and none linear
> equation but in my case none of the equations are linear in a1 or a5.
>
>> I know nothing about the internals of eliminate(), but Maxima is
>> open-source, so one could find out exactly what it does.
>>
> I was just hoping you are someone would know how it got solutions that
> are different from mine.
I'm pretty sure there is a large set (probably infinite) of linear
combinations of polynomials that will solve this. A good and proper way
of solving problems like this is by using Grobner bases. However, I was
unable to get solutions that way -- ran into memory limits on my 2 Gb
machine.
For quartics, there are non-iterative solvers. Google for >>>solve
roots quartic<<< and pick the Wikipedia entry for details. There's
probably some pretty well-optimized code out there somewhere, too, since
this is a common problem that needs a fast solution. Or you can use the
Maxima results of
grind(optimize(solve(c4*x^4+c3*x^3+c2*x^2+c1*x+c0,x)));
(but check it for numerical problems before relying on it). Pretty yucky.
>
> Thanks again
Robert H. Lewis
> Actually, I was going to test Robert H Lewis's
> solution next.
> His solution is the only one that doesn't require
> finding roots.
> The problem is that I am not familiar with Fermat and
> Rober H Lewis
> didn't
> post the script like you did. Also, I liked the idea
> of solving for
> the
> acceleratons first and then substituting. If Fermat
> can calculate any
> of the accelerations without finding roots then I may
> have a answer.
It is quite easy to find resultants for a1 or a5. If you want to pursue this, you will find my email address on my web site.
Robert H. Lewis
Fordham University.
http://home.bway.net/lewis/
pnachtwey
interactive mode which has made experimenting slow. I can't up arrow
to get to previous lines to edit them. Fermat may be good but you
should consider putting a nice user interface on it and charge $600
instead of $60. Students have time, I don't.
Now, if you had posted the script you used to solve the problem so I
could cut and paste like I could with Bob Walton's examples then you
would have sold a copy of Fermat the first day. $60 is nothing.
I am an engineer. Not a mathematician. I need solutions I can use and
solve real time. Like I said above I can solve this problem by other
means. Since a motion controller basically does iterations and
calculates the target position, velocity and acceleration millisecond
by millisecond I can find the acceleration at the sample time just
before a1 is reached. The equations for determining if the target
generator must start ramp to ramp down to x7 ( command position ) is
easy. This is never exact but I can insert a small constant velocity
segment or modify the jerk to get the desired results. This is how I
solved the problem. I am surprised that something that looks so
simple at first turns out to be so complicated.
Peter Nachtwey
pnachtwey
me thinking a little differently about how to solve these kinds of
problems and understanding the problem Mathcad has with solving these
problems. I was too trusting that I could just put in a system of
equations and have a solve function really solve the problem
efficiently. I need another math package. I was impressed with what
you can do with Maxima and it is free. I am surprised no one using
Mathematica tried to solve this problem.
I have found a simpler faster brute force way to solve this problem
Peter Nachtwey
Bob Walton
> ...
> Does any body have a explanation for what the eliminate function
> really does?
Just because I'm interested in this problem:
Maxima's eliminate() function performs "successive resultants". The
resultant() function forms the resultant of two polynomials (see below
for explanation). Example:
q1:2*x^2+y*x+z;
q2:3*x+5*y-a-1;
q3:z^2+x-y^2+5;
eliminate([q1,q2,q3],[y,z]);
gives:
[7425*x^8-1170*x^7+1299*x^6+12076*x^5+22887*x^4-5154*x^3-
1291*x^2+7688*x+15376]
This may also be done "manually" using resultant():
q4:resultant(q3,q2,y);
q5:resultant(q2,q1,y);
q6:resultant(q4,q5,z);
This generates the same result as above (minus the square brackets).
Then the natural question is "OK, so what is a resultant?". The
resultant is the determinant of the Sylvester matrix of the two
polynomials. OK, so what is a "Sylvester matrix?". Given two
polynomials of degree m and n, respectively, with the first polynomial
having coefficients p[m],p[m-1],p[m-2],...,p[1],p[0] and the second
having coefficients q[n],q[n-1],q[n-2],...,q[1],q[0], then the Sylvester
matrix is an (m+n,m+n) square matrix with the coefficients of the two
polynomials arranged as follows (generally, the polynomial coefficients
are zero-padded and shifted right on each successive row:
[p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
[0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0] <--shifted right
[0,0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
...
[0,0,...,0,p[m],p[m-1],p[m-2],...,p[1],p[0]] <--row n
[0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0] <--row n+1
[0,0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
...
[0,0,...,0,q[n],q[n-1],q[n-2],...,q[1],q[0]] <--row m+n
The determinant, naturally, is taken symbolically. Example using q3 and
q2 eliminating y from above:
resultant(q3,q2,y);
gives 24*z^2+(6*x-2)*z-9*x^2+31*x+124
The Sylvester matrix is:
m:matrix([-1,0,z^2+x+5],[5,-z+3*x-1,0],[0,5,-z+3*x-1]);
and the determinant is the same as above (once it is ratsimp()'d).
So that is what eliminate() does (or the equivalent). Note carefully
that changing the order of the equations supplied to eliminate() or the
order of the variables eliminated may change the results, perhaps
drastically. And the solution obtained from eliminate() doesn't
necessarily solve the original equations (in the example case above, it
doesn't). Using a Grobner basis package works better for that:
load(grobner);
ans:poly_grobner([q3,q2,q1],[y,z,x]);
will generate a set of seven equations. (Note: all variables present
in the equation must be listed in the second argument in the order
desired.) The seventh output equation is a quartic in x, and may be
solved for x. The sixth equation is linear in z and does not contain y,
so z may be solved from x. The second equation (which is q2) is linear
in y and contains x and z, and may be used to compute y now that we have
x and z. The resulting set of x, y and z values solves the equations
q1, q2, and q3. One can usually apply solve() directly to the output of
poly_grobner -- in this example, one can, obtaining numeric results.
I note the Maxima version 14 seems to have an improved Grobner basis
package -- it solves your original problem, whereas version 13 would
never terminate.
Peter Nachtwey
"Bob Walton" <s...@sig.invalid> wrote in message
news:480430c0$0$5713$4c36...@roadrunner.com...
> pnachtwey wrote:
>> ...
Hi, I am back after doing taxes :( I got more ram 4GB for my computer so I
can better solve these problems and it was bad and my computer crashed :(
It was a bad week.
>
> Just because I'm interested in this problem:
>
I have a few more problems like this. One problem looks even simpler. It
has only 4 equations with 4 unknowns. Most of the others Mathcad solved
easily but about 1/10 of these result in wrong answers and crashes. I have
a iterative solution for the original problem Hopefully I will be able to
get suitable solutions with what I have learned here.
The whole point of this exercise is that some motion controllers cannot
change the motion profile on-the-fly or have limitations on when the
commands can be issued. Also, some do not generate optimal motion profiles
that execute the command in the minimal time. The natural thing to do is to
determine when the jerk ( third derivative of position ) needs to be 0, or
plus or minus the jerk limit. This results in near sigthed planning and the
results are not optimal. What I have done if given the target generator a
"look ahead" without the recursive search that a chess type program would
do. I have a C program that where I specify the beginning and ending
position velocity and acceleration and velocity, acceleration, deceleration
and jerk limits. The C program generates a series of third order equations
or motion segments that are executed one after the other until the final
position, velocity and acceleration is reached. In some cases geting to the
final position, velocity and acceleration requires 13 segments when the
target generation is going one direction and a command is given to go to a
position behind it. Finding the optimal path required solving the system of
motion segments all at once. The C program graphs the motion profile so I
can see if the profile is right.
With what I have learned I will try doing this with 7th order equations.
This will result in fewer segments but I must be sure the motion segment
doesn't exceed the limits the user specifies.
> Maxima's eliminate() function performs "successive resultants". The
> resultant() function forms the resultant of two polynomials (see below for
> explanation). Example:
>
> q1:2*x^2+y*x+z;
> q2:3*x+5*y-a-1;
> q3:z^2+x-y^2+5;
> eliminate([q1,q2,q3],[y,z]);
>
> gives:
>
> [7425*x^8-1170*x^7+1299*x^6+12076*x^5+22887*x^4-5154*x^3-
> 1291*x^2+7688*x+15376]
I cut and pasted your example above and I found it really should be
q1:2*x^2+y*x+z;
q2:3*x+5*y-z-1;
q3:z^2+x-y^2+5;
eliminate([q3,q2,q1],[y,z]);
This results in the same answer you posted above. What is interesting is
that the answer changes depending on the order of q1,q2 and q3.
This is curious.
eliminate([q1,q2,q3],[y,z]);
[x^2*(45*x^4+3*x^3+11*x^2+81*x+124)]
eliminate([q2,q1,q3],[y,z]);
[25*(45*x^4+3*x^3+11*x^2+81*x+124)]
The order of x is not even the same!!!!!
I will investigate this
>
> This may also be done "manually" using resultant():
>
> q4:resultant(q3,q2,y);
> q5:resultant(q2,q1,y);
> q6:resultant(q4,q5,z);
>
> This generates the same result as above (minus the square brackets).
>
I got
(%i37) q4:resultant(q3,q2,y);
q5:resultant(q2,q1,y);
q6:resultant(q4,q5,z);
(%o37) 24*z^2+(6*x-2)*z-9*x^2+31*x+124
(%o38) (x+5)*z+7*x^2+x
(%o39) 25*(45*x^4+3*x^3+11*x^2+81*x+124)
which is the same result as my example above with [q2,q1,q3]
I modified your solving sequence and got
q4:resultant(q3,q2,z);
q5:resultant(q2,q1,z);
q6:resultant(q4,q5,y);
(%o40) 24*y^2+(30*x-10)*y+9*x^2-5*x+6
(%o41) -(x+5)*y-2*x^2-3*x+1
(%o42) 45*x^4+3*x^3+11*x^2+81*x+124
> Then the natural question is "OK, so what is a resultant?". The resultant
> is the determinant of the Sylvester matrix of the two polynomials. OK, so
> what is a "Sylvester matrix?".
You anticipatted my question
> Given two polynomials of degree m and n, respectively, with the first
> polynomial having coefficients p[m],p[m-1],p[m-2],...,p[1],p[0] and the
> second having coefficients q[n],q[n-1],q[n-2],...,q[1],q[0], then the
> Sylvester matrix is an (m+n,m+n) square matrix with the coefficients of
> the two polynomials arranged as follows (generally, the polynomial
> coefficients are zero-padded and shifted right on each successive row:
>
> [p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
> [0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0] <--shifted right
> [0,0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
> ...
> [0,0,...,0,p[m],p[m-1],p[m-2],...,p[1],p[0]] <--row n
> [0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0] <--row n+1
> [0,0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
> ...
> [0,0,...,0,q[n],q[n-1],q[n-2],...,q[1],q[0]] <--row m+n
>
> The determinant, naturally, is taken symbolically. Example using q3 and
> q2 eliminating y from above:
>
> resultant(q3,q2,y);
>
> gives 24*z^2+(6*x-2)*z-9*x^2+31*x+124
>
> The Sylvester matrix is:
>
> m:matrix([-1,0,z^2+x+5],[5,-z+3*x-1,0],[0,5,-z+3*x-1]);
>
I will have to study this
> and the determinant is the same as above (once it is ratsimp()'d).
>
> So that is what eliminate() does (or the equivalent). Note carefully that
> changing the order of the equations supplied to eliminate() or the order
> of the variables eliminated may change the results, perhaps drastically.
Yes, I can see.
And the solution obtained from eliminate() doesn't
> necessarily solve the original equations (in the example case above, it
> doesn't). Using a Grobner basis package works better for that:
>
> load(grobner);
> ans:poly_grobner([q3,q2,q1],[y,z,x]);
>
> will generate a set of seven equations. (Note: all variables present in
> the equation must be listed in the second argument in the order desired.)
> The seventh output equation is a quartic in x, and may be solved for x.
> The sixth equation is linear in z and does not contain y, so z may be
> solved from x. The second equation (which is q2) is linear in y and
> contains x and z, and may be used to compute y now that we have x and z.
> The resulting set of x, y and z values solves the equations q1, q2, and
> q3. One can usually apply solve() directly to the output of
> poly_grobner -- in this example, one can, obtaining numeric results.
>
Did it! I understand how I can get the right answer but I can't solve for
this real time. However, I can use the answer in my test program to compare
with my real time solutions..
> I note the Maxima version 14 seems to have an improved Grobner basis
> package -- it solves your original problem, whereas version 13 would never
> terminate.
Yes, I downloaded the latest wxMaxima. I will tell others but Maxima. It
will take a while before I know Maxima as well as I know Mathcad..
Bob, I really appreciate the time you have taken to go through these
examples. Maxima is a great teaching tool. I like the way I can just cut
and paste the text and run it and then try different solve orders easily.
Peter Nachtwey
Bob Walton
> "Bob Walton" <s...@sig.invalid> wrote in message
> news:480430c0$0$5713$4c36...@roadrunner.com...
>> pnachtwey wrote:
>>> ...
...
> I have a few more problems like this. One problem looks even simpler. It
> has only 4 equations with 4 unknowns. Most of the others Mathcad solved
> easily but about 1/10 of these result in wrong answers and crashes. I have
> a iterative solution for the original problem Hopefully I will be able to
> get suitable solutions with what I have learned here.
>
...
As I mentioned, I'm also interested in these sorts of problems. I
decided to see if I could come up with a continuous-jerk solution for
the problem starting at an arbitrary acceleration, velocity and position
and ending at a specified position with zero acceleration and velocity.
I assumed a jerk function with enough "humps" to give enough
parameters to solve the constraint equations (final accel=0, final
velocity=0, final position difference from start specified, and the time
of the final stop specified). Without loss of generality, I set the
starting position to zero at time zero. The jerk function I picked is
jerk=j*t*(t1-t)*(t2-t)*(t3-t)
where j, t1 and t2 are constants, and t3 is the specified stopping time.
The jerk is zero at the start and finish. This Maxima code will
solve it:
jerk:j*t*(t1-t)*(t2-t)*(t3-t);
acc:integrate(jerk,t)+a0;
vel:integrate(acc,t)+v0;
pos:integrate(vel,t); /*position starts at 0*/
q1:ratsimp(pos-x3),t=t3; /*position at x3 at time t3*/
q2:ratsimp(vel),t=t3; /*velocity zero at time t3*/
q3:ratsimp(acc),t=t3; /*acceleration zero at time t3*/
ans:solve([q1,q2,q3],[j,t1,t2]);
ans=ans[1];
list:[x3=10,t3=9,a0=-100,v0=62];
ans,list,numer;
Note that two solutions are given by solve(), in which t1 and t2 are
simply interchanged. It doesn't matter which is which. Also, sometimes
t1 and/or t2 may not be in the range of 0 to t3, which is no problem
(not all the "humps" were required). And "j" will be negative if the
"humps" need to go the other way. I think this result is sort of
elegant. Is is possible, though, to choose a bad input relationship
x3=t3*(5*v0+a0*t3)/10, which will cause a divide by zero -- example is
above but with t3=2 and v0=50.
Also, I note that if a smooth (rather than merely continuous) jerk
function is desired,
jerk=j*t^2*(t1-t)*(t2-t)*(t3-t)^2
will give a very similar result (amazingly, only slightly longer
expressions).
********************
Also, I thought about the real-time solution of a quartic a bit more.
If the quartic is a4*x^4+a3*x^3+a2*x^2+a1*x+a0, one may divide through
by a0 (if a0 is 0, then one root is zero, and a cubic remains) and get:
b4*x^4+b3*x^3+b2*x^2+b1*x+1, where b4=a4/a0, etc. Then substitute
y*b4^(-1/4) for x, and get y^4+c3*y^3+c2*y^2+c1*y+1, where
c3=b3/b4^(3/4), etc. This leaves only three coefficients. The desired
root may then be looked up in a three-dimensional linear interpolation
lookup table over the range of coefficient values of interest for a
given problem, rather than be computed by formula (of course, not
forgetting to convert "y" back to "x"). Using a 10x10x10 table should
take less than 10K of ROM, and the computation should be quite fast,
particularly if the intervals are equal permitting indexed lookups.
pnachtwey
> As I mentioned, I'm also interested in these sorts of problems. I
> decided to see if I could come up with a continuous-jerk solution for
> the problem starting at an arbitrary acceleration, velocity and position
> and ending at a specified position with zero acceleration and velocity.
You need to use higher order polynomials. To specify the PVAJ
( positoin, velocity and jerk ) at the two end points requires solving
for 4 equations and four unkowns. Start with the general equation for
a 7th order polynomial. Since you know the first for coefficents
because they are the inital position, velocity, acceleration=2*c and
jerk=6*d.
y(t)=a+b*t+c*t^2+d*t^3+e*t^4+f*t^5+g*t^5+h*t^7
Now differentiate it 4 times and you have four equations and four
unknowns (e,f,g,h). What you don't have this the time interval
between the two points. If you know the time between the two points
then the problem is easy but that is the problem. What is the time
interval between the initial and final point where the velocity,
acceleraton, deceleration and jerk limits are not exceed. A method of
choosing the minimum time that doesn't exceed the user specified
limits is required. This can be done easily enough in Mathcad or
Maxima but not real time.
The order of the polynomial must be equal to the number of specified
items minus one. This means that if you want to specify the beginning
and ending points position, velocity and acceleration, 6 items, then a
5th order polynimial is required.
> I assumed a jerk function with enough "humps" to give enough
> parameters to solve the constraint equations (final accel=0, final
> velocity=0, final position difference from start specified, and the time
> of the final stop specified). Without loss of generality, I set the
> starting position to zero at time zero. The jerk function I picked is
>
> jerk=j*t*(t1-t)*(t2-t)*(t3-t)
???? It looks like you are trying to find a spline that goes between
four points and then integrate the spline over the four points. The
problem with what you are attempting to do is that when the time moves
past time t1 you will them be using t1,t2,t3,and t4 and there will be
a discontinuity in the jerk. The discontinuity usually will not be
too bad but it will be there. What you need to look at is cubic and
quintic splines. Here is a link where I was asking about the fast way
to find the coefficients for quintic splines.
Splines provide an easy way to go make a smooth motion profile between
many points however, you can see that splines dodge the big issue
too. The intervals between points are given. The trick is to
calculate the intervals between the points that are as short as
possible without exceeding the user velocity, acceleration and jerk
limits.
Look at the link to the movies where we execute cubic splines as a
function of a feed chaing positon, velocity and acceleration.
>
> where j, t1 and t2 are constants, and t3 is the specified stopping time.
> The jerk is zero at the start and finish. This Maxima code will
> solve it:
>
> jerk:j*t*(t1-t)*(t2-t)*(t3-t);
> acc:integrate(jerk,t)+a0;
> vel:integrate(acc,t)+v0;
> pos:integrate(vel,t); /*position starts at 0*/
> q1:ratsimp(pos-x3),t=t3; /*position at x3 at time t3*/
> q2:ratsimp(vel),t=t3; /*velocity zero at time t3*/
> q3:ratsimp(acc),t=t3; /*acceleration zero at time t3*/
> ans:solve([q1,q2,q3],[j,t1,t2]);
> ans=ans[1];
> list:[x3=10,t3=9,a0=-100,v0=62];
> ans,list,numer;
I will look at your example. A plot would be good showing the
position, velocity, acceleration and jerk as a function of time would
be nice. Look at this related .pdf generated by a Mathcad worksheet.
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20Cubic%20Interpolation%20NG.pdf
This worksheet shows cubic interpolation, not splines, like you are
doing and the effect of the discontinuity AND resolution of the
points. In the real world there isn't infinite resolution. I like to
use sinewaves because for interpolation problems because:
1. Sinewaves are higher order functions.
2. It is easy to modify the amplitude and frequency of the sinewave
that will challenge the interpolation technique.
>
> Note that two solutions are given by solve(), in which t1 and t2 are
> simply interchanged. It doesn't matter which is which. Also, sometimes
> t1 and/or t2 may not be in the range of 0 to t3, which is no problem
> (not all the "humps" were required). And "j" will be negative if the
> "humps" need to go the other way. I think this result is sort of
> elegant. Is is possible, though, to choose a bad input relationship
> x3=t3*(5*v0+a0*t3)/10, which will cause a divide by zero -- example is
> above but with t3=2 and v0=50.
>
> Also, I note that if a smooth (rather than merely continuous) jerk
> function is desired,
>
> jerk=j*t^2*(t1-t)*(t2-t)*(t3-t)^2
>
> will give a very similar result (amazingly, only slightly longer
> expressions).
>
> ********************
>
> Also, I thought about the real-time solution of a quartic a bit more.
That will not do. You need a higher order equation.
> If the quartic is a4*x^4+a3*x^3+a2*x^2+a1*x+a0, one may divide through
> by a0 (if a0 is 0, then one root is zero, and a cubic remains) and get:
> b4*x^4+b3*x^3+b2*x^2+b1*x+1, where b4=a4/a0, etc. Then substitute
> y*b4^(-1/4) for x, and get y^4+c3*y^3+c2*y^2+c1*y+1, where
> c3=b3/b4^(3/4), etc.
Again I will have to study your technique and see how it can be
applied to higher order polynomials.
> This leaves only three coefficients. The desired
> root may then be looked up in a three-dimensional linear interpolation
> lookup table over the range of coefficient values of interest for a
> given problem, rather than be computed by formula (of course, not
> forgetting to convert "y" back to "x"). Using a 10x10x10 table should
> take less than 10K of ROM, and the computation should be quite fast,
> particularly if the intervals are equal permitting indexed lookups.
The problem is that a the order of the polynomial is usually higher
and one needs to find the roots of higher order polynomials.
Meanwhile I am still digesting what has been presented above.
Peter Nachtwey
Bob Walton
Correction to the Sylvester matrix description (correction to row n+1
and thereafter):
> pnachtwey wrote:
>> ...
...
> polynomials. OK, so what is a "Sylvester matrix?". Given two
> polynomials of degree m and n, respectively, with the first polynomial
> having coefficients p[m],p[m-1],p[m-2],...,p[1],p[0] and the second
> having coefficients q[n],q[n-1],q[n-2],...,q[1],q[0], then the Sylvester
> matrix is an (m+n,m+n) square matrix with the coefficients of the two
> polynomials arranged as follows (generally, the polynomial coefficients
> are zero-padded and shifted right on each successive row:
>
> [p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
> [0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0] <--shifted right
> [0,0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
> ...
> [0,0,...,0,p[m],p[m-1],p[m-2],...,p[1],p[0]] <--row n
> [0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0] <--row n+1
> [0,0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
> ...
> [0,0,...,0,q[n],q[n-1],q[n-2],...,q[1],q[0]] <--row m+n
>
[p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
[0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0] <--shifted right
[0,0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
...
[0,0,...,0,p[m],p[m-1],p[m-2],...,p[1],p[0]] <--row n
[q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0] <--row n+1
[0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
[0,0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
...
[0,0,...,0,q[n],q[n-1],q[n-2],...,q[1],q[0]] <--row m+n
...
Bob Walton
Correction to the Sylvester matrix description (correction to row n+1
and thereafter):
> pnachtwey wrote:
>> ...
...
> polynomials. OK, so what is a "Sylvester matrix?". Given two
> polynomials of degree m and n, respectively, with the first polynomial
> having coefficients p[m],p[m-1],p[m-2],...,p[1],p[0] and the second
> having coefficients q[n],q[n-1],q[n-2],...,q[1],q[0], then the Sylvester
> matrix is an (m+n,m+n) square matrix with the coefficients of the two
> polynomials arranged as follows (generally, the polynomial coefficients
> are zero-padded and shifted right on each successive row:
>
> [p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
> [0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0] <--shifted right
> [0,0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
> ...
> [0,0,...,0,p[m],p[m-1],p[m-2],...,p[1],p[0]] <--row n
> [0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0] <--row n+1
> [0,0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
> ...
> [0,0,...,0,q[n],q[n-1],q[n-2],...,q[1],q[0]] <--row m+n
>
[p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
[0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0] <--shifted right
[0,0,p[m],p[m-1],p[m-2],...,p[1],p[0],0,0,...,0]
...
[0,0,...,0,p[m],p[m-1],p[m-2],...,p[1],p[0]] <--row n
[q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0] <--row n+1
[0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
[0,0,q[n],q[n-1],q[n-2],...,q[1],q[0],0,0,...,0]
...
[0,0,...,0,q[n],q[n-1],q[n-2],...,q[1],q[0]] <--row m+n
...
pnachtwey
Sylverter Matrix.
http://planetmath.org/encyclopedia/DerivationOfSylvestersMatrixForTheResultant.html
This makes sense. I also noticed the term resultant. I haven't had
time to persue how this resultant relates to the Dixon Resultant that
Rober H Lewis mentioned.
Peter Nachtwey