Pagina 2 van 2

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 11 Jan 2018, 20:11
door Simotion
Arie,

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

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 13 Jan 2018, 23:07
door Simotion
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?

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 16 Jan 2018, 10:18
door arie
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?

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 16 Jan 2018, 20:16
door Simotion
Arie,

Bedankt, ik bekijk dit verder.

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 21 Jan 2018, 16:28
door Simotion
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.

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 23 Jan 2018, 13:07
door arie
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...

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 23 Jan 2018, 18:34
door Simotion
Arie,

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

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 01 Apr 2018, 20:31
door Simotion
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

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 08 Apr 2018, 09:12
door arie
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.

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 09 Apr 2018, 22:48
door Simotion
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?

Re: Coëfficienten van polynoom berekenen

BerichtGeplaatst: 10 Apr 2018, 17:25
door arie
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]));
}