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: 16
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: 16
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: 2989
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: 16
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: 16
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: 2989
Geregistreerd: 09 Mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

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

Re: Coëfficienten van polynoom berekenen

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

Re: Coëfficienten van polynoom berekenen

Berichtdoor 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.
arie
Moderator
Moderator
 
Berichten: 2989
Geregistreerd: 09 Mei 2008, 09:19

Re: Coëfficienten van polynoom berekenen

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

Re: Coëfficienten van polynoom berekenen

Berichtdoor 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: Alles selecteren
\\ 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]));
}
arie
Moderator
Moderator
 
Berichten: 2989
Geregistreerd: 09 Mei 2008, 09:19

Vorige

Terug naar Praktijkproblemen

Wie is er online?

Gebruikers in dit forum: Geen geregistreerde gebruikers en 2 gasten

cron

Wie is er online?

Er zijn in totaal 2 gebruikers online :: 0 geregistreerd, 0 verborgen en 2 gasten (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 2 gasten
Copyright © 2009 Afterburner - Free GPL Template. All Rights Reserved.