Numerik: Octave

21/02/2011 - 20:46 von 0liver ojo Bedf0rd | Report spam
Hi!

Ich versuche gerade meine ersten Schritte in Octave. Und zwar möchte ich
den pH-Wert einer Lösung einer Sàure und einer Base berechnen. In
etwa so:

function y = f(x)
KA = 10^-4.75;
KB1 = 10^-4.07;
KB2 = 10^-7.15;
KW = 10^-14;
HA0 = 0.002;
B0 = 0.001;
y = zeros(6, 1);
y(1) = x(2) + x(6) - x(5) - x(4);
# a oh h bh
y(2) = KA - x(2)*x(5)/x(1);
# a h ha
y(3) = KB1 - x(4)*x(6)/x(3);
# bh oh b
y(4) = KW - x(5)*x(6);
# h oh
y(5) = HA0 - x(1) - x(2);
# ha a
y(6) = B0 - x(3) - x(4);
# b bh
endfunction;

Aufrufe von fsolve scheitern mit verschiedenen info-Werten und
zum Teil physikalisch unsinnigen Lösungen (Konzentrationen < 0,
h * oh != 10^-14 = Ionenprodukt von Wasser).

Wie macht man es richtig?

Ciao und Danke,
Oliver
 

Lesen sie die antworten

#1 Christian Gollwitzer
21/02/2011 - 22:45 | Warnen spam
Am 21.02.11 20:46, schrieb 0liver 'ojo' Bedf0rd:
Ich versuche gerade meine ersten Schritte in Octave. Und zwar möchte ich
den pH-Wert einer Lösung einer Sàure und einer Base berechnen. In
etwa so:

function y = f(x)
KA = 10^-4.75;
KB1 = 10^-4.07;
KB2 = 10^-7.15;
KW = 10^-14;
HA0 = 0.002;
B0 = 0.001;
y = zeros(6, 1);
y(1) = x(2) + x(6) - x(5) - x(4);
# a oh h bh
y(2) = KA - x(2)*x(5)/x(1);
# a h ha
y(3) = KB1 - x(4)*x(6)/x(3);
# bh oh b
y(4) = KW - x(5)*x(6);
# h oh
y(5) = HA0 - x(1) - x(2);
# ha a
y(6) = B0 - x(3) - x(4);
# b bh
endfunction;

Aufrufe von fsolve scheitern mit verschiedenen info-Werten und
zum Teil physikalisch unsinnigen Lösungen (Konzentrationen< 0,
h * oh != 10^-14 = Ionenprodukt von Wasser).



Ich kenne diese Gleichungen nicht, daher kann ich nur allgemeine Tipps
geben. Leider löst die Numerik meist nur Probleme, deren Lösung man
schon etwas kennt. fsolve verwendet ein iteratives Verfahren um sich der
Lösung anzunàhern. Dabei spielt der Startpunkt eine wesentliche Rolle -
wenn Du zu weit von der Lösung entfernt bist, konvergiert der
Algorithmus nicht zu einer sinnvollen Lösung.

Bei einem kurzen Test ist mir aufgefallen, dass Dein Problem
augenscheinlich ganz miserabel skaliert ist:

res=fsolve('f', [0.1 0.2 0.5 0.5 0.5 0.5 ]'/10)
f(res)


ans
-2.6196e-13
8.3536e-13
3.4327e-13
-5.1561e-14
2.3980e-13
2.2598e-13


Für fsolve stellt das eine Lösung dar, weil der Betrag von res sehr
klein ist. Deine Gleichung 4 (h*oh-10^-14=0) ist mit einem Residuum von
-5.156*10^-14 offensichtlich schlecht erfüllt, wàhrend das für eine
Konzentrationssumme eine prima Genauigkeit wàre. Die Lösung heißt
skalieren und umformulieren. Multipliziere Gleichung 4 einfach mit 10^10
durch um auf dieselbe Skala zu kommen wie Gleichung 6 etc.
Möglicherweise hilft hier auch Logarithmieren, da die Skalen so
unterschiedlich sind. Also z.B. Gleichung 4 :

y(4)=log(KW)-(log(x(5)*x(6)))


Weiterhin hilft es, die Anzahl der Dimensionen des Problems zu
reduzieren. Wenn Du eine der ganz simplem Gleichungen (z.B. 4, 5 oder 6)
auflöst und in die anderen einsetzt, kannst Du eine Unbekannte
eliminieren. Je weniger Unbekannte das System hat, umso weniger leicht
kann sich der Algorithmus "verirren". Allerdings musst Du hier auch
etwas aufpassen, wenn z.B. ha und a sehr unterschiedliche
Größenordnungen haben, so kann (wegen der endlichen Genauigkeit) HA0-ha
sehr unterschiedlich zu a sein (s. Fließkommaarithmetik). In diesem Fall
müsste man so auflösen, dass de kleinere der beiden Zahlen als Variable
frei bleibt.

So nun wünsch ich viel Spaß beim Spielen ;)

Christian

Ähnliche fragen