4 onbekenden, 4 vergelijkingen en veel machtsverheffen

[Marcel L.]

Voor wiskundige talenten wellicht een aardigheidje, maar voor mij een onoverkoombare horde: 4 onbekenden bepalen wanneer je 4 vergelijkingen hebt.
In de praktijk gaat het om het bepalen van de calibratieconstanten van een meetinstrument. Op 4 punten (n=2, n=6, n=10 en n=14) geeft het instrument de ruwe waarde Xn. Met behulp van de onderstaande formule is het mogelijk de lijn die de Xn waarden verbindt om te vormen tot een mooie rechte lijn binnen de meettoleranties.

De onbekenden zijn: A, B, C, en V. Bekend zijn n en Xn. Voorbeeld waarden van X2 = 1,33; X6 = 12,13; X10 = 23,56 en X14 = 32,16
De vergelijking is n*C+V+(B*A^Xn) = Xn

En als dit echt te makkelijk was; nog tips om dit in Excel te automatiseren?

[op=op]

Deze gelijkheden van elkaar aftrekken geeft

Evenzo is
Deze gelijkheden van elkaar aftrekken geeft

Evenzo is
Deze gelijkheden op elkaar delen geeft voor n=2

Deze transcendente gelijkheid is op te lossen naar A met b.v. Newton-Raphson of de halveringsmethode.
Als je A gevonden hebt is de rest eenvoudig.
Met Excel bemoei ik me niet.

[Marcel L.]

Hulde,op=op!

Met mijn roestige kennis van 23 jaar geleden kan ik dit allemaal (nog net) volgen en met behulp van Wikipedia valt de halveringsmethode goed in een Excel-macro te vangen.

Bedankt!

[SafeX]

Evenzo is
Deze gelijkheden van elkaar aftrekken geeft

[op=op]

Goed opgemerkt SaveX.
Fouten maken blijft onmenselijk.
Dus het moet zijn

[arie]

Deze transcendente gelijkheid is op te lossen naar A met b.v. Newton-Raphson of de halveringsmethode.

Niet voor A>0, terwijl ik vermoed dat je juist A>0 wilt: het gaat om calibratie van een meetinstrument dat wsch reele waarden afgeeft.

Ga uit van

Dan moet gelden A ongelijk nul EN A ongelijk 1, anders is links de noemer nul.
Dit is overigens ook direct te zien in het stelsel: de term (B*A^x) reduceert dan naar nul resp B voor alle x>0.

Uitgewerkt vind je nu met gegeven waarden:
0.63*A^32.16 + 1.57*A^23.56 - 5.03*A^12.13 + 2.83*A^1.33 = 0
ofwel (deel door A^1.33, ga na dat dit mag):
0.63*A^30.83 + 1.57*A^22.23 - 5.03*A^10.8 + 2.83 = 0
definieer de functie
g(A) = 0.63*A^30.83 + 1.57*A^22.23 - 5.03*A^10.8 + 2.83
dan zie je direct
g(1) = 0
maar A moest ongelijk 1 zijn.

Bepaal nu g'(A) = de afgeleide van g(A):
g'(A) = A^9.8 * (19.4229*A^20.03 + 34.9011*A^11.43 - 54.324)
dan zie je dat g'(1)=0.
Omdat A>0 is de eerste factor altijd positief, en is
g'(A) negatief voor 0<A<1, dus g(A) dalend op dit gebied
g'(A) positief voor A>1, dus g(A) stijgend op dit gebied,
waardoor er geen ander nulpunt dan A=1 is voor g(A), en A=1 was niet toegestaan...

Het gegeven stelsel heeft dus geen oplossingen.

[Marcel L.]

Arie - tot die conclusie was ik inmiddels ook gekomen. Door (bij gebrek aan kunde) te itereren, ben ik er na heel veel tijd achter gekomen dat:
bereik A: 0,85 tot 0,98
bereik B: -4 tot -9
bereik C: 1,5 tot 2,5
bereik V: 6,0 tot 10,0

Om d.m.v. trial & error de kromme grafiek van Xn dusdanig te bewerken m.b.v. de formule tot er een mooie rechte lijn uitkomt kost me teveel tijd en daarom ben ik opzoek naar een manier om in iedergeval 1 van de variabelen A, B, C (of V) te berekenen, waardoor het calibreren van het meetinstrument sneller uitgevoerd kan worden.

Dus opnieuw de vraag: Wie voelt zich uitgedaagd?

[arie]

Je zou hier met de brutekracht-methode de optimale waarden voor A, B, C en V kunnen zoeken.

Echter: waarom gebruik je hier die ingewikkelde vergelijking: n*C+V+(B*A^Xn) = Xn ?

Als ik je punten plot liggen deze vrijwel op een lijn, en kan je ook heel goed uitkomen met lineaire regressie. Losjes gezegd: in het algemeen wordt daarbij de lijn gezocht die zo goed mogelijk door je puntenverzameling loopt.
Zie bijvoorbeeld http://phoenix.phys.clemson.edu/tutoria ... ssion.html
Het gaat in jouw geval om de lijn L bepaald door a en b:
L: Xn = a*n + b
In Excel zijn hiervoor standaardfuncties:
a = SLOPE ~= 2.598
b = INTERCEPT ~= -3.489
En de correlatie coefficient (in Excel de functie CORREL) van je punten ~= 0.998329486, wat behoorlijk hoog is (1=perfecte correlatie).
Je lijn zou daarmee dus worden:
Xn = 2.598*n - 3.489
Ben je wellicht naar zoiets op zoek, of moeten we toch uitgaan van je oorspronkelijke vergelijking?

[Marcel L.]

Waarom die ingewikkelde formule? Zo gaat dat met productontwikkeling - afhankelijk van de gebruikte onderdelen ziet de 'ruwe' curve er soms heel anders uit. Over het algemeen in het lage bereik (n=1,5 en n=2) een grotere afwijking naar beneden; in het midden mooi lineair en vervolgens in het hogere bereik (n=10 t/m n=14) weer een grotere afwijking waarbij ook naar boven moet worden gecorrigeerd.

De gebruikte formule is door derden bedacht en geeft voldoende vrijheid om de ruwe data om te vormen tot het gewenste resultaat. Deze formule wordt ook gebruikt in de on-board software van het meetinstrument, alleen wordt dáár de A, B, C en V gebruikt die ik handmatig bepaald heb.

Inmiddels heb ik tientallen keren de variabelen iteratief bepaald (samen met een grafiekje voor de visuele feedback) maar dat is natuurlijk niet dé manier om zo'n probleem op te lossen.

[op=op]

Merk op dat
,
waaruit volgt dat .

[arie]

Dan moeten we het dus op de formule n*C+V+(B*A^Xn) = Xn houden.
Zou je s.v.p. eerst nog wel even kunnen nakijken of dit mogelijk
n*C+V+(B*A^n) = Xn,
zou moeten zijn, dus links alleen termen met 'n' erin (A^n i.p.v. A^Xn), en rechts de 'Xn'.

[Marcel L.]

De formule n*C+V+(B*A^Xn) = Xn is correct en juist die A^Xn is essentieel omdat je hiermee vooral het lage bereik (n=<4) hier goed mee kunt beïnvloeden.

Situatieschets: n is de waarde die ik wil bereiken (is dus steeds een bekende) en Xn is de ruwe meetwaarde die door de formule omgevormd moet worden tot n. Dus ik stel mijn mechanische opstelling dusdanig in dat n=2 en doe dan metingen X2. Vervolgens pas ik de opstelling aan naar n =4 en vervolgens meet ik X4 enzovoort.

[arie]

Je zou dan kunnen itereren over A: van 0.85 tot 0.98 met stapgrootte bv 0.0001.
Voor elke A krijg je dan een stelsel van 4 vergelijkingen in de vorm n*C+V+(A^Xn)*B = Xn met 3 onbekenden (C, V en B; voor elke A heb je immers steeds 4 waarden van n, Xn en A^Xn).
Pas hierop een lineair regressiemodel toe, en je vindt de waarden van C, V en B die het beste passen gegeven de waarde van A.
Selecteer dan die A waarbij het model het allerbeste past (= zo groot mogelijke correlatiecoefficient R),
en je hebt je antwoord.

De formules vind je bv terug op http://www.ncssm.edu/courses/math/Stat_ ... ession.pdf

Als ik dit vertaal in Pari/GP (http://pari.math.u-bordeaux.fr/), dan krijg ik:

Code: Selecteer alles

\\main:
{
x=matrix(4,3);
y=matrix(4,1);

y=[1.33; 12.13; 23.56; 32.16];
x=[2,1,0; 6,1,0; 10,1,0; 14,1,0];

\\bepaal l=aantal meetwaarden (in dit voorbeeld: l=4):
l=length(y[,1]);

\\bereken y average:
yave=0;
for(i=1,l, yave+=y[i,1]; );
yave/=l;

\\bereken noemer R^2 = (correlatie coefficient)^2:
y0=0.0;
for(i=1,l, y0+=(y[i,1]-yave)^2; );

\\iterate:
r2max=0;
forstep(a=0.85, 0.98, 0.0001,
  \\bepaal 3e kolom van x, afhankelijk van de huidige a:
  for(i=1,l,x[i,3]=a^y[i,1]);
  \\los regressiemodel op:
  r=(x~ * x)^(-1)*(x~ * y);
  \\neem oplossingen over (voor de duidelijkheid):
  c=r[1,1];
  v=r[2,1];
  b=r[3,1];
  \\bereken teller R^2:
  y1=0.0;
  for(i=1,l, y1+=(x[i,1]*c + v + x[i,3]*b - yave)^2; );
  r2=y1/y0;
  if(r2>r2max, r2max=r2; rmax=r; amax=a; ); \\bewaar beste resultaat
  );

\\toon resultaat:
print("R^2 = ",r2max);
print("C = ",rmax[1,1]);
print("V = ",rmax[2,1]);
print("B = ",rmax[3,1]);
print("A = ",amax);
}

In dit voorbeeld krijg ik als uitkomst:
R^2 = 0.999357369939
C = 1.61256957468
V = 23.2516562527
B = -26.0133785365
A = 0.980000000000
waarbij A dus zo ver mogelijk naar de 1 getrokken wordt.

Excel kent ook matrixbewerkingen (MINVERSE, MMULT, etc), dus dit zou ook daarin te implementeren moeten zijn.

[Marcel L.]

Dit ziet er goed uit; het code-voorbeeld is ook bijzonder paraktisch. Ik nog wel even tijd nodig om het allemaal te doorgronden en praktisch uit te werken maar ik ga ervoor... Alvast bedankt Arie!

[op=op]

Als je mijn bijdragen eens wat beter leest scheelt dat een hoop werk. Succes!