Coëfficienten van polynoom berekenen

Wiskunde is niet alleen een vak op school. Kom je ergens in de praktijk (bijvoorbeeld tijdens je werk) een wiskundig probleem tegen dan kun je hier om hulp vragen.
Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 11 jan 2018, 20:11

Arie,

Alvast hartelijk bedankt om dit uit te werken.
Ik bekijk hoe ik hiermee verder aan de slag kan.

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 13 jan 2018, 23:07

Arie,

Uiteindelijk zou ik tot de oplossing moeten komen waarbij de coëfficienten A tem F gevonden worden aan de hand van de doorgegeven variabelen (Xts,Vts,Ats,Xte,Vte,Ate).
De motion controller moet dus een functie hebben die op elke moment met nieuwe variabelen de coëfficienten kan berekenen.
Hoe vorm je de vermenigvuldiging van de inverse matrix met de andere matrix om tot formules voor de coëfficienten op basis van de ingegeven parameters?

arie
Moderator
Moderator
Berichten: 3071
Lid geworden op: 09 mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

Bericht door arie » 16 jan 2018, 10:18

Als je software niet met matrices en vectoren kan rekenen, maar alleen met basale rekenkundige bewerkingen (zoals optellen, aftrekken, vermenigvuldigen en delen van getallen), dan wordt je programma-code wat langer.

Een mogelijke oplossing:
Werk eerst via Gauss-eliminatie achtereenvolgens F, E en D weg, zie https://nl.wikipedia.org/wiki/Gauss-eliminatie en met name de op die wiki-pagina genoemde elementaire rijoperaties die we kunnen gebruiken:
- twee rijen (=vergelijkingen) verwisselen
- een rij met een scalair (=getal) ongelijk aan 0 vermenigvuldigen of erdoor delen
- bij een rij (een veelvoud van) een andere rij optellen of aftrekken.
Dit levert:

Code: Selecteer alles

Ga uit van het eerder beschreven stelsel van 6 vergelijkingen met 6 onbekenden:

[1]  xs =    A*s5 +    B*s4 +   C*s3 +   D*s2 + E*s + F
[2]  vs =  5*A*s4 +  4*B*s3 + 3*C*s2 + 2*D*s  + E
[3]  as = 20*A*s3 + 12*B*s2 + 6*C*s  + 2*D
[4]  xe =    A*e5 +    B*e4 +   C*e3 +   D*e2 + E*e + F
[5]  ve =  5*A*e4 +  4*B*e3 + 3*C*e2 + 2*D*e  + E
[6]  ae = 20*A*e3 + 12*B*e2 + 6*C*e  + 2*D

waarbij s5=ts^5, s4=ts^4 etc. 


Trek vergelijking [1] van vergelijking [3] af, dat geeft vergelijking [7]:

[2]  vs =       5*s4*A +    4*s3*B +    3*s2*C +     2*s*D +       E
[3]  as =      20*s3*A +   12*s2*B +     6*s*C +       2*D
[7]  xe-xs = (e5-s5)*A + (e4-s4)*B + (e3-s3)*C + (e2-s2)*D + (e-s)*E
[5]  ve =       5*e4*A +    4*e3*B +    3*e2*C +     2*e*D +       E
[6]  ae =      20*e3*A +   12*e2*B +     6*e*C +       2*D


Deel [7] door (e-s), dat geeft [8]:

[2]  vs =                              5*s4*A +              4*s3*B +        3*s2*C +   2*s*D + E
[3]  as =                             20*s3*A +             12*s2*B +         6*s*C +     2*D
[8] (xe-xs)/(e-s) = (e4+e3*s+e2*s2+e*s3+s4)*A + (e3+e2*s+e*s2+s3)*B + (e2+e*s+s2)*C + (e+s)*D + E
[5]  ve =                              5*e4*A +              4*e3*B +        3*e2*C +   2*e*D + E
[6]  ae =                             20*e3*A +             12*e2*B +         6*e*C +     2*D


Trek [2] af van [8] en van [5], dat geeft resp. [9] en [10]:

[3]  as =                                  20*s3*A +               12*s2*B +           6*s*C +       2*D
[9]  (xe-xs)/(e-s)-vs= (e4+e3*s+e2*s2+e*s3-4*s4)*A + (e3+e2*s+e*s2-3*s3)*B + (e2+e*s-2*s2)*C +   (e-s)*D
[10] ve-vs =                           5*(e4-s4)*A +           4*(e3-s3)*B +     3*(e2-s2)*C + 2*(e-s)*D 
[6]  ae =                                  20*e3*A +               12*e2*B +           6*e*C +       2*D


Deel [3] en [6] door 2, dat geeft [11] en [12]:

[11] as/2 =                                10*s3*A +                6*s2*B +           3*s*C +         D
[9]  (xe-xs)/(e-s)-vs= (e4+e3*s+e2*s2+e*s3-4*s4)*A + (e3+e2*s+e*s2-3*s3)*B + (e2+e*s-2*s2)*C +   (e-s)*D
[10] ve-vs =                           5*(e4-s4)*A +           4*(e3-s3)*B +     3*(e2-s2)*C + 2*(e-s)*D 
[12] ae/2 =                                10*e3*A +                6*e2*B +           3*e*C +         D


Trek [11] af van [12], dat geeft [13]:

[9]  (xe-xs)/(e-s)-vs= (e4+e3*s+e2*s2+e*s3-4*s4)*A + (e3+e2*s+e*s2-3*s3)*B + (e2+e*s-2*s2)*C +   (e-s)*D
[10] ve-vs =                           5*(e4-s4)*A +           4*(e3-s3)*B +     3*(e2-s2)*C + 2*(e-s)*D 
[13] ae/2-as/2 =                      10*(e3-s3)*A +           6*(e2-s2)*B +       3*(e-s)*C


Vermenigvuldig [11] met (e-s), dat levert [14]:
[14] (e-s)*as/2 =                    10*(e-s)*s3*A +          6*s2*(e-s)*B +     3*s*(e-s)*C +   (e-s)*D
ofwel:
[14] (e-s)*as/2 =                (10*e*s3-10*s4)*A +       (6*e*s2-6*s3)*B + (3*e*s-3*s^2)*C +   (e-s)*D


Trek [14] af van [9], dat levert [15]:

[15] (xe-xs)/(e-s)-vs-(e-s)*as/2 = (e4+e3*s+e2*s2-9*e*s3+6*s4)*A + (e3+e2*s-5*e*s2+3*s3)*B + (e2-2*e*s+s2)*C
[10] ve-vs =                                       (5*e4-5*s4)*A +           (4*e3-4*s3)*B +   (3*e2-3*s2)*C + 2*(e-s)*D 
[13] ae/2-as/2 =                                    10*(e3-s3)*A +             6*(e2-s2)*B +       3*(e-s)*C


Vermenigvuldig [14] met 2, dat levert [16]:

[16] (e-s)*as =                                (20*e*s3-20*s4)*A +       (12*e*s2-12*s3)*B + (6*e*s-6*s^2)*C +   2*(e-s)*D


Trek [16] af van [10], dat levert [17]:

[15] (xe-xs)/(e-s)-vs-(e-s)*as/2 = (e4+e3*s+e2*s2-9*e*s3+6*s4)*A + (e3+e2*s-5*e*s2+3*s3)*B +     (e2-2*e*s+s2)*C
[17] ve-vs-(e-s)*as =                     (5*e4-20*e*s3+15*s4)*A +   (4*e3-12*e*s2+8*s3)*B + (3*e2-6*e*s+3*s2)*C
[13] ae/2-as/2 =                                    10*(e3-s3)*A +             6*(e2-s2)*B +           3*(e-s)*C

We houden nu een stelsel over van 3 vergelijkingen met 3 onbekenden (A, B en C) in de vorm

p = m1*A + m2*B + m3*C
q = m4*A + m5*B + m6*C
r = m7*A + m8*B + m9*C

waarbij p, q, r en m1 t/m m9 bekende constanten zijn.
We kunnen dit stelsel oplossen door verder te gaan met de Gauss eliminatie, maar dat levert nogal wat werk. Als alternatief kunnen we het vanaf dit punt ook direct oplossen via de regel van Cramer (zie https://nl.wikipedia.org/wiki/Regel_van_Cramer):

Herschrijf het stelsel in matrixvorm (v = M * x):



en bereken de 4 benodigde determinanten via het schema van Sarrus
(https://nl.wikipedia.org/wiki/Regel_van_Sarrus).

In onderstaand programma heb ik bovenstaande verwerkt in de functie solvems().
(alles wat op een regel achter \\ staat is commentaar, dus geen code)
Je zal dit waarschijnlijk wel kunnen omzetten naar jouw programmeertaal.

Code: Selecteer alles

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ solve master slave probleem:                        \\
\\ s = starttijd                                       \\
\\ e = eindtijd                                        \\
\\ x, v en a zijn positie, snelheid en versnelling     \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
solvems(s,e,xs,vs,as,xe,ve,ae)={
\\ bereken eerst de benodigde machten:
s2=s*s;
s3=s*s2;
s4=s*s3;
s5=s*s4;
e2=e*e;
e3=e*e2;
e4=e*e3;
e5=e*e4;

\\ bereken matrix M (zie vergelijkingen [15], [17] en [13]):
m1=e4+e3*s+e2*s2-9*e*s3+6*s4;
m2=e3+e2*s-5*e*s2+3*s3;
m3=e2-2*e*s+s2;

m4=5*e4-20*e*s3+15*s4;
m5=4*e3-12*e*s2+8*s3;
m6=3*e2-6*e*s+3*s2;

m7=10.0*(e3-s3);
m8=6*(e2-s2);
m9=3*(e-s);

\\ bereken vector v (zie eveneens vergelijkingen [15], [17] en [13]):
p=(xe-xs)/(e-s)-vs-(e-s)*as/2;
q=ve-vs-(e-s)*as;
r=ae/2.0-as/2;

\\ bereken de determinant van M:
detm=m1*m5*m9+m2*m6*m7+m3*m4*m8-m7*m5*m3-m8*m6*m1-m9*m4*m2;

\\ gebruik de regel van Cramer om A, B en C op te lossen:
A=(p*m5*m9+m2*m6*r+m3*q*m8-r*m5*m3-m8*m6*p-m9*q*m2)/detm;
B=(m1*q*m9+p*m6*m7+m3*m4*r-m7*q*m3-r*m6*m1-m9*m4*p)/detm;
C=(m1*m5*r+m2*q*m7+p*m4*m8-m7*m5*p-m8*q*m1-r*m4*m2)/detm;

\\ gebruik vergelijking [11] om D te vinden:
D = as/2 - 10*s3*A - 6*s2*B - 3*s*C;
\\ gebruik vergelijking [2] om E te vinden:
E = vs - 5*s4*A - 4*s3*B - 3*s2*C - 2*s*D;
\\ gebruik vergelijking [1] om F te vinden:
F = xs - A*s5 - B*s4 - C*s3 - D*s2 - E*s;

\\ output: print de resultaten:
print("A = ", A);
print("B = ", B);
print("C = ", C);
print("D = ", D);
print("E = ", E);
print("F = ", F);

}

\\ MAIN PROGRAM:
\\ roep solvems() aan met de gegevens van voorbeeld 2:
{
solvems(700,1000,600,1,0,1000,1,0);
}
Kom je hiermee verder?

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 16 jan 2018, 20:16

Arie,

Bedankt, ik bekijk dit verder.

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 21 jan 2018, 16:28

Bovenstaande heb ik geïmplementeerd in de motion controller, en dit werkt zoals het moet.
Ik ben nu bezig met een uitbreiding waarbij de eindpositie van de master niet gekend is.
De eindpositie van de slave is de eindpositie van de master plus een gekende offset.
Wel gekend is de maximale acceleratie van de slave (tov de master) overheen het volledige traject.
De maximale acceleratie zal optreden op de masterpositie waar de eerste afgeleide van de acceleratie van de slave 0 is.
Dit is waar positie master = (-24*B + SQRT(576B²-1440AC))/120A of positie master = (-24*B - SQRT(576B²-1440AC))/120A.
Daarmee worden de vergelijkingen wel een stuk moeilijker.

arie
Moderator
Moderator
Berichten: 3071
Lid geworden op: 09 mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

Bericht door arie » 23 jan 2018, 13:07

Dit wordt inderdaad lastig.
We hadden:

Code: Selecteer alles

[1]  xs =    A*s5 +    B*s4 +   C*s3 +   D*s2 + E*s + F
[2]  vs =  5*A*s4 +  4*B*s3 + 3*C*s2 + 2*D*s  + E
[3]  as = 20*A*s3 + 12*B*s2 + 6*C*s  + 2*D
[4]  xe =    A*e5 +    B*e4 +   C*e3 +   D*e2 + E*e + F
[5]  ve =  5*A*e4 +  4*B*e3 + 3*C*e2 + 2*D*e  + E
[6]  ae = 20*A*e3 + 12*B*e2 + 6*C*e  + 2*D
Nu wordt
e = s + dt (waarbij dt de tijdsduur van de verplaatsing is)
xe = eindpositie master + offset = s + dt + offset
waarbij dt ook een onbekende is (in het stelsel ontstaan dan producten van onbekenden).
Bovendien moeten we de maximale acceleratie van de slave begrenzen.

Je zou dit numeriek kunnen oplossen.

We hadden al de functie
solvems(s,e,xs,vs,as,xe,ve,ae)
waarmee we de polynoomconstanten A t/m F kunnen bepalen voor gegeven e en xe.
Daarmee kunnen we de maximale acceleratie voor die gegeven e en xe berekenen, zie code hieronder:
aslave(xm) bepaalt de acceleratie als functie van de master positie
getamax() berekent de maximale acceleratie, met als input gegeven:
- dt = tijdsduur verplaatsing
- offset = offset eindpositie slave t.o.v. master
- s, xs, vs, as, ve, ae als in solvems()
De formules voor pm1 en pm2 had je al bepaald.
De versnellingen neem ik hieronder absoluut: de vertraging = negatieve versnelling verwacht ik eveneens begrensd (zo niet, pas dit dan aan).

Code: Selecteer alles

aslave(xm)={20*A*xm^3 + 12*B*xm^2 + 6*C*xm + 2*D}

getamax(s,dt,xs,vs,as,offset,ve,ae)={
\\ calculate A..F:
solvems(s,s+dt,xs,vs,as,s+dt+offset,ve,ae);
\\ find extreme acceleration master positions:
pm1=(-24*B-sqrt(576*B^2-1440*A*C))/(120*A);
pm2=(-24*B+sqrt(576*B^2-1440*A*C))/(120*A);
\\ calculate absolute value of accelerations:
a1 = abs(aslave(pm1));
a2 = abs(aslave(pm2));
\\ get maximum value of those two:
amax=max(a1,a2);
\\return value:
amax
}
In het eerdere voorbeeld 2 is dt=300 en offset=0.
We krijgen dan:
pm1 = 763.397459621...
pm2 = 936.602540378...
amax = a1 = a2 = 0.006415002990...

We hebben nu dus de functie getamax() waar we dt in kunnen stoppen en amax uit krijgen.
Via het bisectie algoritme (https://en.wikipedia.org/wiki/Bisection_method) kunnen we nu dt vinden voor een gegeven amax.
Dit gebeurt in de nieuwe functie finddt() hieronder:

Code: Selecteer alles

finddt(atarget,s,xs,vs,as,offset,ve,ae)={
\\ determine start interval:
dthigh=100;
ahigh=getamax(s,dthigh,xs,vs,as,offset,ve,ae);
while(ahigh<atarget,
  dthigh=dthigh/2;
  ahigh=getamax(s,dthigh,xs,vs,as,offset,ve,ae);
  );
while(ahigh>atarget,
  dtlow=dthigh;
  dthigh=dthigh*2;
  ahigh=getamax(s,dthigh,xs,vs,as,offset,ve,ae);
  );
\\ now:  dtlow <= dttarget <= dthigh

\\ iterations bisection method:
for(i=1,50,
  dtmid=(dtlow+dthigh)/2.0;
  amid=getamax(s,dtmid,xs,vs,as,offset,ve,ae);
  if(amid<atarget,
    dthigh = dtmid
    ,\\ELSE:
    dtlow = dtmid
    );
  );
\\return value:
(dtlow+dthigh)/2.0
}
Eerst bepalen we het interval [dtlow, dthigh] waarbinnen de tijd met de gewenste amax moet liggen. Daarna halveren we dit interval 50 keer: for(i=1,50,...), zodat het oorspronkelijke interval ongeveer een factor 10^15 kleiner wordt (elke 10 halveringen leveren een factor 2^10 ~= 1000 kleiner).

Stel in voorbeeld 2 wil je de maximale versnelling brengen naar
amax = 0.01,
dan vind je via
dt = finddt(0.01,700,600,1,0,0,1,0)
een waarde van
dt = 240.2811414134753853...
en hierbij hoort een
amax = 0.01000000000000000355864879...

Ter controle:
solvems(700, 700+240.2811414134753853, 600, 1, 0, 700+240.2811414134753853+0, 1, 0)
ofwel
solvems(700, 940.2811414134753853, 600, 1, 0, 940.2811414134753853, 1, 0)
levert
A = 0.0000000007491224610518077656
B = -0.000003071928613681327500339
C = 0.0050027844015779016163304406
D = -4.043867160498191179397450613
E = 1623.6854978560383829807966872
F = -258774.9414804240429206434824
waarbij
pm1=(-24*B-sqrt(576*B^2-1440*A*C))/(120*A) = 750.77737986860741811...
en
aslave(pm1) = 0.01000000000000000356
aslave(pm1-0.1) = 0.0099999688081030...
aslave(pm1+0.1) = 0.0099999688380679...
dus de slave heeft een maximum acceleratie bij pm1=750.777 van ongeveer 0.0100000...

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 23 jan 2018, 18:34

Arie,

Alvast hartelijk bedankt.
Zodra ik de gelegenheid zie probeer ik dit te implementeren in de motion controller.

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 01 apr 2018, 20:31

De bovenstaande logica (met opgegeven max acceleratie en te zoeken masterafstand) heb ik succesvol geïmplementeerd.
Ik ben nu ondertussen bezig met het uitrekenen van een polynoom die de positie van de slave voorstelt in functie van de positie van de master zoals in de post van 16 Jan 2018, 10:18, maar nu echter met extra voorwaarden startpunt en eindpunt.
In het voorbeeld van 16 jan was voor startpunt en eindpunt van de curve telkens de positie, snelheid en acceleratie opgegeven. Nu zou ook de eerste afgeleide van de acceleratie voor startpunt en eindpunt opgegeven zijn : dit noemt men de 'jerk' (j) of wijziging van de acceleratie.
Men heeft dan 8 vergelijkingen voor 8 onbekenden.
Om dit tot code om te vormen die door mijn controller kan worden uitgevoerd (zoals in voorbeeld van 16 jan) wordt dit echter een heel stuk complexer. Na een aantal pogingen loop ik dan ook vast in de veelheid
van vergelijkingen. Hoe kan men deze vergelijkingen het beste herleiden om er de coëfficiënten (A tem H) uit te berekenen?
[1]  xs =    A*s7 +    B*s6 +   C*s5 +   D*s4 + E*s3 + F*s2 + Gs + H
[2]  vs =  7*A*s6 +  6*B*s5 + 5*C*s4 + 4*D*s3  + 3*E*s2 + 2*F*s + G
[3]  as = 42*A*s5 + 30*B*s4 + 20*C*s3 + 12*D*s2 + 6*E*s + 2*F
[4]  js = 210*A*s4 + 120*B*s3 + 60*C*s2 + 24*D*s + 6*E
[5]  xe =    A*e7 +    B*e6 +   C*e5 +   D*e4 + E*e3 + F*e2 + Ge + H
[6]  ve =  7*A*e6 +  6*B*e5 + 5*C*e4 + 4*D*e3  + 3*E*e2 + 2*F*e + G
[7]  ae = 42*A*e5 + 30*B*e4 + 20*C*e3 + 12*D*e2 + 6*E*e + 2*F
[8]  je = 210*A*e4 + 120*B*e3 + 60*C*e2 + 24*D*e + 6*E

arie
Moderator
Moderator
Berichten: 3071
Lid geworden op: 09 mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

Bericht door arie » 08 apr 2018, 09:12

In een programmeertaal zonder matrices en vectoren wordt dat inderdaad een monnikenwerk.
Werk eerst de variabelen weg waarvan er de minste zijn, te beginnen met H:

[1] xs = A*s7 + B*s6 + C*s5 + D*s4 + E*s3 + F*s2 + Gs + H
[2] vs = 7*A*s6 + 6*B*s5 + 5*C*s4 + 4*D*s3 + 3*E*s2 + 2*F*s + G
[3] as = 42*A*s5 + 30*B*s4 + 20*C*s3 + 12*D*s2 + 6*E*s + 2*F
[4] js = 210*A*s4 + 120*B*s3 + 60*C*s2 + 24*D*s + 6*E
[5] xe = A*e7 + B*e6 + C*e5 + D*e4 + E*e3 + F*e2 + Ge + H
[6] ve = 7*A*e6 + 6*B*e5 + 5*C*e4 + 4*D*e3 + 3*E*e2 + 2*F*e + G
[7] ae = 42*A*e5 + 30*B*e4 + 20*C*e3 + 12*D*e2 + 6*E*e + 2*F
[8] je = 210*A*e4 + 120*B*e3 + 60*C*e2 + 24*D*e + 6*E

gaat over in:

[5-1] xe-xs = A*(e7-s7) + B*(e6-s6) + C*(e5-s5) + D*(e4-s4) + E*(e3-s3) + F*(e2-s2) + G*(e-s)
[2] vs = 7*A*s6 + 6*B*s5 + 5*C*s4 + 4*D*s3 + 3*E*s2 + 2*F*s + G
[3] as = 42*A*s5 + 30*B*s4 + 20*C*s3 + 12*D*s2 + 6*E*s + 2*F
[4] js = 210*A*s4 + 120*B*s3 + 60*C*s2 + 24*D*s + 6*E
[6] ve = 7*A*e6 + 6*B*e5 + 5*C*e4 + 4*D*e3 + 3*E*e2 + 2*F*e + G
[7] ae = 42*A*e5 + 30*B*e4 + 20*C*e3 + 12*D*e2 + 6*E*e + 2*F
[8] je = 210*A*e4 + 120*B*e3 + 60*C*e2 + 24*D*e + 6*E

een stelsel van 7 vergelijkingen met 7 onbekenden.

Bepaal vervolgens om G te elimineren:
[6] - [2]
en
[5-1] - (e-s)*[2]
Daarna vervalt [2] en hou je 6 vergelijkingen met 6 onbekenden over.


PS:
Voordat je dit doet zou ik wel eerst controleren of je getallen niet te groot worden voor de nauwkeurigheid van je computer.
Bijvoorbeeld: als s=100, dan is s^7 = 10^14, een getal van 15 cijfers.

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 09 apr 2018, 22:48

Los van de uitwerking van de uitwerking van de formules voor toepassing in mijn controller,
welke online matrix calculatoren zouden er aan te bevelen zijn om reeds voor specifieke data,
de matrix voor de zevende orde polynoom uit te rekenen?

arie
Moderator
Moderator
Berichten: 3071
Lid geworden op: 09 mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

Bericht door arie » 10 apr 2018, 17:25

We hebben:



of kortweg:



vermenigvuldig links en rechts met de inverse matrix:



ofwel, de oplossing van het stelsel is:



Je kan deze oplossing uit laten rekenen door Pari/GP, de webversie vind je hier:
http://pari.math.u-bordeaux.fr/gp.html

Copy/paste onderstaande code naar het blauwe invoerveld van die pagina, en klik vervolgens op de knop "Evaluate with PARI" (die knop staat direct onder dat invoerveld).
In het roze veld krijg je dan opnieuw de code te zien, en het resultaat van de berekening.

In de code hieronder is alles wat volgt op \\ commentaar.
Je kan alle getallen zelf naar wens aanpassen.

Code: Selecteer alles

\\ reken met 60 cijfers nauwkeurigheid:
\p 60


\\ MAIN program:
{

\\ startwaarden:
xs=12;
vs=0.4;
as=0.02;
js=0.005;

\\ eindwaarden:
xe=36;
ve=0.3;
ae=0.03;
je=0.003;

\\ starttijd en eindtijd:
s=2.3;
e=5.6;

\\ vector v:
v=[xs,vs,as,js, xe,ve,ae,je];

\\ matrix M:
M=[ s^7,     s^6,     s^5,    s^4,    s^3,   s^2, s, 1;
    7*s^6,   6*s^5,   5*s^4,  4*s^3,  3*s^2, 2*s, 1, 0;
    42*s^5,  30*s^4,  20*s^3, 12*s^2, 6*s,   2,   0, 0;
    210*s^4, 120*s^3, 60*s^2, 24*s,   6,     0,   0, 0;
    e^7,     e^6,     e^5,    e^4,    e^3,   e^2, e, 1;
    7*e^6,   6*e^5,   5*e^4,  4*e^3,  3*e^2, 2*e, 1, 0;
    42*e^5,  30*e^4,  20*e^3, 12*e^2, 6*e,   2,   0, 0;
    210*e^4, 120*e^3, 60*e^2, 24*e,   6,     0,   0, 0 ];

\\ bereken oplossing:
w = M^-1 * v~;

\\ print oplossing in nette vorm naar beeldscherm:
print("A = ", w[1]);
print("B = ", w[2]);
print("C = ", w[3]);
print("D = ", w[4]);
print("E = ", w[5]);
print("F = ", w[6]);
print("G = ", w[7]);
print("H = ", w[8]);
print();

\\ ter controle:
\\ M*w moet weer de oorspronkelijke vector v opleveren:
c=M*w;
\\ print controlevector c:
for(i=1,8,print(c[i]));
}

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 10 mei 2018, 21:48

In mijn voorgaande vraag had ik het over zevende orde polynomen.
Het specifieke nut van een zevende orde polynoom is dat je als je dit gebruikt om een positie van een slave tov een master te definiëren, je voor begin en eindsituatie zowel positive, snelheid, acceleratie als jerk (afgeleide van acceleratie) kan vooropstellen.
Nu zijn er controllers waar ik mee werk, die maar polynomen tot de zesde orde kunnen verwerken.
De fabricant heeft echter wel een functieblok waarmee je at runtime in de controller een polynoomcurve kan berekenen waaraan je voor begin- en eindsituatie de positive, snelheid, acceleratie en jerk kan meegeven.Ik zie dat de functie in feite twee polynomen oplevert van de zesde orde, die dan na elkaar gezet worden. Op die manier worden de voorwaarden voor begin- en eindsituatie toch voldaan.
Dit is echter een afgeschermde functieblok waarin ik de code niet kan zien. Daarom had ik graag gekeken of ik dezelfde functionaliteit ook zelf zou kunnen programmeren.
Ik weet echter totaal niet waar ik het tussenpunt (eindpunt eerste polynoom = beginpunt tweede polynoom) zou leggen en met welke criteria voor de verschillende condities (snelheid, acceleratie, etc.).

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 26 mei 2018, 22:13

Betreffende mijn voorgaande post om een polynoomcurve te maken waarbij er gekenden zijn voor start en eindpunt, betreffende de positie, snelheid, acceleratie en afgeleide van de acceleratie (jerk), maar waar dit moet uitgevoerd worden met polynoom maximum tot de zesde orde.
Ik heb een reeks testen kunnen doen met de functieblok van de fabrikant. Hieruit kan ik de coëfficiënten van de polynoom halen. Ik zag ook dat wat normaal één polynoom zou zijn, is opgesplitst in twee trajecten met elk een andere polynoom. Uit de testen bleek de eerste polynoom tot de helft van de masterpositie te lopen en de tweede vanaf de helft van de masterpositie tot de eindpositie.
Voorbeeld1 :
positie van de master aan begin = px0 = 0 (mm)
positie van de slave aan begin = py0 = 0 (mm)
positie van de master aan eind = pxt = 250 (mm)
positie van de slave aan eind = pyt = 150 (mm)

snelheid van master aan begin = vx0 = 100 mm/sec
acceleratie van master aan begin = ax0 = 0
snelheid van slave aan begin = vy0 = 0 (factor tov master : 0 = stilstand slave)
acceleratie van slave aan begin = ay0 = 0
jerk van slave aan begin = jy0 = 0.0

snelheid van master aan eind = vxt = 100 mm/sec
acceleratie van master aan eind = axt = 0
snelheid van slave aan eind = vyt = 1 (factor tov master : 1 = zelfde snelheid als master)
acceleratie van slave aan eind = ayt = 0
jerk van slave aan eind = jyt = 0.0

Resultaat : (coëfficiënten voor bijhorende machten : [1] xs = A*s6 + B*s5 + C*s4 + D*s3 + E*s2 + F*s + G
curve 1: masterpositie 0-125
A: 0 , B: -1.0, C: 2.0, D: 0, E: 0, F: 0, G: 0
curve 2: masterpositie 125-250
A: 0 , B: -0.052631578, C: 0.263157894, D: -0.526315789, E: 0.526315789, F: 0.789473684, G: 2.775557561
Voorbeeld2 :
positie van de master aan begin = px0 = 110 (mm)
positie van de slave aan begin = py0 = 100 (mm)
positie van de master aan eind = pxt = 250 (mm)
positie van de slave aan eind = pyt = 350 (mm)

snelheid van master aan begin = vx0 = 100 mm/sec
acceleratie van master aan begin = ax0 = 0
snelheid van slave aan begin = vy0 = 1.0 (factor tov master : 1 = zelfde snelheid als master)
acceleratie van slave aan begin = ay0 = 0.1
jerk van slave aan begin = jy0 = 0.0

snelheid van master aan eind = vxt = 100 mm/sec
acceleratie van master aan eind = axt = 0
snelheid van slave aan eind = vyt = 1.5 (factor tov master : 1.5 = 50% sneller dan master)
acceleratie van slave aan eind = ayt = 0.4
jerk van slave aan eind = jyt = 0.0

Resultaat : (coëfficiënten voor bijhorende machten : [1] xs = A*s6 + B*s5 + C*s4 + D*s3 + E*s2 + F*s + G
curve 1: masterpositie 110-180
A: 0 , B: -0.766247379, C: 0.885744234, D: 0.0, E: 0.684835779, F: 0.195667365, G; -5.551115123
curve 2: masterpositie 180-250
A: 0 , B: -4.384144360, C: 13.09380410, D: -8.533772826, E: -3.445356968, F: 2.645541957, G: 0.640285927

Behalve de vaststelling dat ze de curves verdelen over telkens de helft van de masterafstand heb ik niet echt een idee welke criteria je zou kunnen nemen om dan met het respecteren van de 8 voorwaarden (start en eindpunt) te komen tot 2 aaneensluitende polynomen, die in deze voorbeelden van de vijfde orde blijken te zijn, maar mag tot zesde orde gaan (niet zevende orde). In het middenpunt waar de eerste curve eindigt en de tweede begint zie ik dat er continuïteit is in positie, snelheid en acceleratie.
Ik had graag geweten hoe ik dit zou aanpakken als ik het zelf programmeer.
Een zevende order polynoom opstellen voor de 8 voorwaarde lukt me nu al, maar deze curve uit twee deelcurves laten bestaan om zo aan de start en eindpunten aan de voorwaarden te voldoen, maar er tevens voor zorgen dat de maximale orde van de polynomen zes is, lukt me nog niet. Advies, tips zijn dus welkom.

Simotion
Nieuw lid
Nieuw lid
Berichten: 19
Lid geworden op: 22 dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Bericht door Simotion » 05 jun 2018, 21:10

Geen ideëen omtrent bovenstaande?

arie
Moderator
Moderator
Berichten: 3071
Lid geworden op: 09 mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

Bericht door arie » 08 jun 2018, 12:45

Hoe zijn positie en tijd nu gedefinieerd???
Wat opvalt is dat je constante coefficienten nu ongeveer dezelfde grootte-orde hebben:
A: 0 , B: -1.0, C: 2.0, D: 0, E: 0, F: 0, G: 0
curve 2: masterpositie 125-250
A: 0 , B: -0.052631578, C: 0.263157894, D: -0.526315789, E: 0.526315789, F: 0.789473684, G: 2.775557561
De term B*x^5 is met deze getallen veruit dominant in het gebied van 50 tot 250

Vergelijk:
In je post van 9 Jan 2018 kwamen je bijvoorbeeld uit op
y = 0.000000000246914 x5 -0.000001049382716 x4 + 0.001765432099 x3 -1.469135802 x2 + 605.9382716 x -98765.43210
(je functie FCCalcCamPoly5Order)
In mijn berekening van 11 Jan 2018 kwam ik uit op dezelfde getallen.


Als ik voor het huidige probleem een soortgelijk stelsel opstel, kom ik voor je voorbeeld 1 uit op:

Eerste helft: uitgaande van x(t) = B*t^5 + C*t^4 + D*t^3 + E*t^2 + F*t + G:

B = 0.0000002625536
C = -0.000050432
D = 0
E = 0
F = 0
G = 0

Hiermee wordt:
x(0) = 0
v(0) = 0
a(0) = 0
j(0) = 0

x(125) = -4300.00
v(125) = -73.5000
a(125) = 0.800000
j(125) = 0.094848


Tweede helft: uitgaande van x(t) = Q*t^5 + R*t^4 + S*t^3 + T*t^2 + U*t + V:

Q = 0.0000003444736
R = -0.00035456
S = 0.139264
T = -25.312
U = 2076.0
V = -64250.0

Nu is:
x(125) = -4300.00
v(125) = -73.5000
a(125) = 0.800000
j(125) = 0.094848
dus precies aansluitend de eindwaarden van de eerste helft

x(250) = 150
v(250) = 100
a(250) = 0
j(250) = 0


Mogelijk werk je nu met andere eenheden voor plaats en tijd.
Of is er iets heel anders aan de hand?

Plaats reactie