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.

Re: Coëfficienten van polynoom berekenen

Berichtdoor 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: 12
Geregistreerd: 22 Dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Berichtdoor 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?
Simotion
Nieuw lid
Nieuw lid
 
Berichten: 12
Geregistreerd: 22 Dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Berichtdoor 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: Alles selecteren
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: Alles selecteren
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ 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?
arie
Moderator
Moderator
 
Berichten: 2964
Geregistreerd: 09 Mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

Berichtdoor Simotion » 16 Jan 2018, 20:16

Arie,

Bedankt, ik bekijk dit verder.
Simotion
Nieuw lid
Nieuw lid
 
Berichten: 12
Geregistreerd: 22 Dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Berichtdoor 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.
Simotion
Nieuw lid
Nieuw lid
 
Berichten: 12
Geregistreerd: 22 Dec 2017, 21:46

Re: Coëfficienten van polynoom berekenen

Berichtdoor arie » 23 Jan 2018, 13:07

Dit wordt inderdaad lastig.
We hadden:

Code: Alles selecteren
[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: Alles selecteren
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: Alles selecteren
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...
arie
Moderator
Moderator
 
Berichten: 2964
Geregistreerd: 09 Mei 2008, 09:19

Vorige

Terug naar Praktijkproblemen

Wie is er online?

Gebruikers in dit forum: Geen geregistreerde gebruikers en 1 gast

Wie is er online?

Er is in totaal 1 gebruiker online :: 0 geregistreerd, 0 verborgen en 1 gast (Gebaseerd op de gebruikers die actief waren gedurende 5 minuten)
De meeste gebruikers ooit tegelijkertijd online was 649 op 31 Okt 2014, 18:45

Gebruikers in dit forum: Geen geregistreerde gebruikers en 1 gast
Copyright © 2009 Afterburner - Free GPL Template. All Rights Reserved.