Solving set of equations doesn't give the correct answer

9 Ansichten (letzte 30 Tage)
Bob photonics
Bob photonics am 23 Mai 2018
Bearbeitet: Bob photonics am 25 Mai 2018
Dear everyone,
I've been trying to translate a set of chemical equations to matlab code, to be able to solve for different chemical species. I've the approximate solution (as it's from a graph) but after entering all the data and checking multiple times I still haven't been able to find what wrong. I'm wondering what is going wrong and if anyone could please help me out. The source for the graph/equation is the article at this link: The chemistry of co-injected BOE The graph I want to reproduce later on is figure 2 in the paper, see the image below:
Now the results I get for 10cc, 40cc and 90cc are respectively: HF 43%, H2F2 48%, F- 3%, HF2- 6% in comparison ~28%, 63%, 2%, 7% (10cc). HF 35%, H2F2 33%, F- 14%, HF2- 18% in comparison ~24%, 44%, 6%, 26% (40cc). HF 21%, H2F2 12%, F- 37%, HF2- 30% in comparison ~18%, 23%, 20%, 45% (90cc).
the script is the following:
clear all;
%Units to be used
%Volume is in CC also cm^3, 1 litre is 1000 CC, 1 cc = 1 ml
%density is in g/cm^3
%weigth percentages are in fractions of 0 to 1
%Molecular weight is in g/mol
% pts=10; %number of points for linear spacing
%weight percentages of NH4OH and HF
dh2o=1.00; %0.997 at 25C when rounded 1
%HF values
dhf=dh2o+(dhf49-dh2o)*xhf/0.49; %@ 25C
%NH4OH (NH3) values
% Vnh3=linspace(0.1*Vhf,1.9*Vhf,pts);
dnh3=0.9; %for ~20-31% @~20-25C
Mnh3=17.03; %The wt% of NH4OH actually refers to the wt% of NH3 dissolved in H2O
if max(nnh3)>=nhf
error(['There are more mols NH4OH,',num2str(max(nnh3)),', than mols HF,',num2str(nhf),'.'])
%%Calculations for species
Vt=(Vhf+Vh2o+Vnh3)/1000; %litre
A=nhf/Vt; %mol/l
B=nnh3/Vt; %mol/l
syms HF F H2F2 HF2 NH3 NH4 H OH
eq2= H*F/HF==6.85*10^(-4);
eq3= NH3*H/NH4==6.31*10^(-10);
eq4= H*OH==10^(-14);
eq5= HF2/(HF*F)==3.963;
eq6= H2F2/(HF^2)==2.7;
eq7= H+NH4==OH+F+HF2;
eq8= HF+F+2*H2F2+2*HF2==A;
eq9= NH3+NH4==B;
varias=[HF, F, H2F2, HF2, NH3, NH4, H, OH];
assume(HF> 0 & F>= 0 & H2F2>= 0 & HF2>= 0& NH3>= 0 & NH4>= 0 & H>= 0 & OH>= 0)
[HF, F, H2F2, HF2, NH3, NH4, H, OH]=vpasolve(eqns,varias);% [0 max([A,B])])
HFf=double(HF)/totalHF %fraction of species for HF
H2F2f=double(H2F2)/totalHF %fraction of species for H2F2
Ff=double(F)/totalHF %fraction of species for F-
HF2f=double(HF2)/totalHF %fraction of species for HF2-
an extra function needed is called mols.m
%%%%amount of mol, Vol=volume, d=density, pwt=%weight, M=molecularweight
function mol=mols(Vol, d, pwt, M)
  2 Kommentare
John D'Errico
John D'Errico am 23 Mai 2018
Bearbeitet: John D'Errico am 23 Mai 2018
It does not give you the answer that you EXPECT. There is no conclusion we can draw so far that it is not correct.
Have you verified that the solution found does NOT solve the equations posed? If it does, that just says there MAY be multiple solutions. Or perhaps it may say that the solution you think is valid, may possibly not be a solution at all.
There is no requirement that a nonlinear system of equations has a unique solution. Or your problem may be that you have formulated the system improperly.
Note that I cannot test your code at all, since I lack at least the function mols.
Bob photonics
Bob photonics am 23 Mai 2018
The article, explains the assumptions made and they say that there can be two solutions, 1 positive and 1 negative but a negative solutions is unrealistic, solve and vpasolve also show all solutions possible. All the variables need to be 0 or greater (except HF has to just be greater). When I run and check it, solve/vpasolve only give 1 solution.
So, yes I might have formulated it inproperly but I've checked this over 5 times by reading it and keeping the article next to it (if you have a better way to check please tell me).
Are there any other ways to check/verify if there are more solutions? Or if the systems is formulated correctly.
Ps this is the equation set:

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Steven Lord
Steven Lord am 23 Mai 2018
Try calling vpasolve with one output argument rather than multiple output arguments. This way you can eliminate the possibility that vpasolve returns the answers in a different order than you expected. When you do this, that output argument will be a struct array you can pass into the subs function directly to test that the values vpasolve returned satisfy your equations. A simpler example of what I mean using solve is:
syms x y
eq1 = x + y == 4;
eq2 = x - y == 2;
S = solve(eq1, eq2)
% Substitute back into one equation
subs(eq1, S) % Displays 4 == 4
% Is that equation always true?
isAlways(subs(eq1, S)) % true
% Are both equations always true?
isAlways(subs([eq1, eq2], S)) % [true true]
% Is the system satisfied?
all(isAlways(subs([eq1, eq2], S))) % Yes!
If you wanted to manually check rather than trusting isAlways:
S.x + S.y % 4
S.x - S.y % 2
  2 Kommentare
Bob photonics
Bob photonics am 23 Mai 2018
Going to try this tomorrow, thanks for the suggestion :D
Bob photonics
Bob photonics am 25 Mai 2018
Bearbeitet: Bob photonics am 25 Mai 2018
So I've found a problem when I try to do it step by step at least.
It's not giving me a solution because matlab won't solve this (2016 both a and b versions)
syms x
>> eq12 = -3 == sqrt(x);
>> solve(eq12)
ans =
Empty sym: 0-by-1
Obviously the answer here would be 9, so what do I change to let matlab solve this and before you tell me to change the equation, this is the equation I would have to change by hand and not just once but a 100 times or so:
93659574124777211691008 H
H + ---------------------------------------------
1208925819614629174706176 H + 762832192176831
== -----------------
100000000000000 H
29250045579840375 #1
- --------------------------------------------------------------------
H (922337203685477580800 H + 927343445063259249) 4611686018427387904
12879770070323045125 #1
+ ----------------------------------------------------------------------
55340232221128654848 H (922337203685477580800 H + 927343445063259249)
/ 2
| 181939 H 852912375078609598437 H
#1 == 9223372036854775808 H - sqrt| --------- + -----------------------
\ 5539 25544128856069301600256
39917248404619332215368770561441 |
+ -------------------------------------- | 9223372036854775808 + 6318009845245521
85070591730234615865843651857942052864 /
The trouble maker is the minus in front of the sqrt I assume
PS when not pretty
H + (93659574124777211691008*H)/(1208925819614629174706176*H + 762832192176831) == 1/(100000000000000*H) - (29250045579840375*(9223372036854775808*H - 9223372036854775808*((181939*H^2)/5539 + (852912375078609598437*H)/25544128856069301600256 + 39917248404619332215368770561441/85070591730234615865843651857942052864)^(1/2) + 6318009845245521))/(4611686018427387904*H*(922337203685477580800*H + 927343445063259249)) + (12879770070323045125*(9223372036854775808*H - 9223372036854775808*((181939*H^2)/5539 + (852912375078609598437*H)/25544128856069301600256 + 39917248404619332215368770561441/85070591730234615865843651857942052864)^(1/2) + 6318009845245521)^2)/(55340232221128654848*H*(922337203685477580800*H + 927343445063259249)^2)

Melden Sie sich an, um zu kommentieren.




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by