2 Gnu Octave | ||
→ | 2.16 Beispiel-Aufgaben | |
→ | 2.16.6 Permanentmagnet |
Dieses Beispiel zeigt die Spline-Interpolation.
Aufgabenstellung |
B-Wert für gegebenes H finden |
H-Wert für vorgegebenes B finden |
Optimalen Arbeitspunkt finden |
Plot |
Gegeben ist die B-H-Kennlinie eines permanentmagnetischen Materials in Form von Messdaten, siehe Tabelle:
H / Acm-1 | B / T |
---|---|
-525 | 0 |
-470 | 0,5 |
-420 | 0,72 |
-330 | 0,94 |
-220 | 1,12 |
-120 | 1,22 |
0 | 1,30 |
Folgende Größen sollen ermittelt werden:
Mit der Octave-Datei
# Vektor mit den H-Werten
h = [ -525 -470 -420 -330 -220 -120 0 ]';
# Vektor mit den B-Werten
b = [ 0 0.50 0.72 0.94 1.12 1.22 1.30 ]';
# Spline-Interpolation
global pp;
pp = interp1(h, b, 'spline', 'pp');
function [y] = Bfct(x)
global pp;
y = ppval(pp, x);
endfunction
# Problem 1: Suche nach B, wenn H=-250 Ampere pro Zentimeter
printf('B1 = %g\n', Bfct(-250));
finden wir die Lösung B1=1,078 T.
Die Funktion Bfct berechnet B für einen gegebenen H-Wert. Dazu wird das vorher mit interp1() angelegte stückweise Polynom pp benutzt.
Gesucht wird die Nullstelle von f(H)=B(H)-0,8 T.
# Vektor mit den H-Werten
h = [ -525 -470 -420 -330 -220 -120 0 ]';
# Vektor mit den B-Werten
b = [ 0 0.50 0.72 0.94 1.12 1.22 1.30 ]';
# Spline-Interpolation
global pp;
pp = interp1(h, b, 'spline', 'pp');
# Erste Ableitung des Splines
global ppd;
ppd = ppder(pp);
function [y, yder] = Bsoll(x)
global pp;
global ppd;
y = ppval(pp, x) - 0.8;
yder = NA;
if (nargout > 1)
yder = ppval(ppd, x);
endif
endfunction
# Problem 2: Suche nach H für vorgegebenes B=0.8 Tesla
[H2, distance, info] = fsolve(@Bsoll, -250, optimset("jacobian", "on"));
if (1 == info)
printf('Loesung: H2 = %g\n', H2);
else
printf('Leider keine Loesung gefunden!\n');
endif
Die Funktion Bsoll implementiert die Funktion, deren
Nullstelle wir suchen (Differenz zum Sollwert). Werden zwei
Argumente in der Rückgabeliste erwartet, wird als zweiter
Rückgabewert die erste Ableitung dieser Funktion zurückgegeben.
Das für die Berechnung der ersten Ableitung verwendete stückweise
Polynom ppd wurde vorher mit der Funktion ppder() aus
dem stückweisen Polynom pp der Funktion gebildet.
Für die Nullstellensuche wird fsolve() benutzt.
Die Funktion f1(H) - im Octave-Code BH genannt - liefert das Produkt B*H. Am Extremwert dieser Funktion ist die erste Ableitung 0. Um die Nullstelle der ersten Ableitung zu suchen, wollen wir Octave auch deren erste Ableitung (also die zweite Ableitung der Originalfunktion) zur Verfügung stellen.
# Vektor mit den H-Werten
h = [ -525 -470 -420 -330 -220 -120 0 ]';
# Vektor mit den B-Werten
b = [ 0 0.50 0.72 0.94 1.12 1.22 1.30 ]';
global pp; # Ergebnis der Spline-Interpolation
pp = interp1(h, b, 'spline', 'pp');
# pp = csape(h, b, 'not-a-knot');
# pp = csape(h, b, 'variational');
global ppd1; # Erste Ableitung von pp
ppd1 = ppder(pp);
global ppd2; # Zweite Ableitung von pp
ppd2 = ppder(ppd1);
# Berechnung von B*H.
function [y] = BH(x)
global pp;
y = x .* ppval(pp, x);
endfunction
# Erste und zweite Ableitung von B*H
function [y, yder] = BHder(x)
global pp; global ppd1; global ppd2;
y = ppval(pp, x) + x .* ppval(ppd1, x);
yder = NA;
if (nargout > 1)
yder = 2 * ppval(ppd1, x) + x * ppval(ppd2, x);
endif
endfunction
# Problem 3: Suche nach Maximum fuer B*H
[H3, diff, info] = fsolve(@BHder, -250, optimset("jacobian", "on"));
if (1 == info)
printf('H3 = %g\nB3 = %g\nBH = %g\n', H3, ppval(pp, H3), BH(H3));
else
printf('Leider keine Loesung fuer Problem 3 gefunden!\n');
end
Im Quelltext implementiert die Funktion BH die Funktion
zur Produktbildung, BHder implementiert die erste und -
falls zwei Rückgabewerte vom Funktionsaufruf erwartet werden - auch
die zweite Ableitung. Mit fsolve() wird nach der Nullstelle
der ersten Ableitung gesucht.
Dabei werden die stückweisen Polynome ppd1 und ppd2
verwendet, die zur Berechnung der ersten und zweiten Ableitung des
stückweisen Polynomes pp dienen.
# Vektor mit den H-Werten
h = [ -525 -470 -420 -330 -220 -120 0 ]';
# Vektor mit den B-Werten
b = [ 0 0.50 0.72 0.94 1.12 1.22 1.30 ]';
global pp; # Ergebnis der Spline-Interpolation
pp = interp1(h, b, 'spline', 'pp');
# pp = csape(h, b, 'not-a-knot');
# pp = csape(h, b, 'variational');
# Berechnung von B
function [y] = B(x)
global pp;
y = ppval(pp, x);
endfunction
# Berechnung von B*H.
function [y] = BH(x)
global pp;
y = abs(x .* ppval(pp, x));
endfunction
graphics_toolkit("gnuplot");
fplot(@B, [-525, 0]);
Um den Plot für |BH| zu erhalten, wird die fplot()-Anweisung ersetzt durch:
fplot(@BH, [-525, 0]);