Solving set of equations doesn't give the correct answer
6 views (last 30 days)
Show older comments
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/190790/image.png)
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:
clc;
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
xhf=0.49;
xnh3=0.28;
%H2O
Vh2o=1800;
dh2o=1.00; %0.997 at 25C when rounded 1
mh2o=18.02;
%HF values
Vhf=100;
dhf49=1.15;
dhf=dh2o+(dhf49-dh2o)*xhf/0.49; %@ 25C
Mhf=20.01;
nhf=mols(Vhf,dhf,xhf,Mhf);
%NH4OH (NH3) values
% Vnh3=linspace(0.1*Vhf,1.9*Vhf,pts);
Vnh3=10;
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
nnh3=mols(Vnh3,dnh3,xnh3,Mnh3);
if max(nnh3)>=nhf
error(['There are more mols NH4OH,',num2str(max(nnh3)),', than mols HF,',num2str(nhf),'.'])
end
n2=nnh3;
n3=nhf-nnh3;
%%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;
eqns=[eq2,eq5,eq6,eq8,eq4,eq3,eq9,eq7];
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])])
totalHF=double(HF)+double(F)+double(H2F2)+double(HF2);
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)
mol=(Vol*d*pwt)/M;
end
2 Comments
John D'Errico
on 23 May 2018
Edited: John D'Errico
on 23 May 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.
Answers (1)
Steven Lord
on 23 May 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 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!