# How to set up and minimize a large system of equations to resolve the individual components

I am trying to set up a system of 400*N equations to deconvolute spectral data for four different compounds.

I monitored 400 different wavelengths in N solvents (3 solvents in this case). In order to do this, I want to try to minimize the following function:

$A_{\lambda_i, N_i} - (\epsilon_{1 (\lambda_i, N_i)}C_1 + \epsilon_{2 (\lambda_i, N_i)}C_2 + \epsilon_{3 (\lambda_i, N_i)}C_3 + \epsilon_{4 (\lambda_i, N_i)}C_4 = 0$

Here $A$ is the absorbance of the spectrum I want to deconvolute at specified wavelength ($\lambda$) in solvent ($N$).
$\epsilon$ is the molar extinction coefficient for each different compound that corresponds to the absorbance at the given wavelength and solvent. Then unknown quantities are $C_1$ to $C_4$

My data is organized as follows.

Each row has 5 entries and I have 400 rows one (for each wavelength) [$A \; \epsilon_1 \; \epsilon_2 \; \epsilon_3 \; \epsilon_4$]

Each solvent has a table organized just like this one. How do I used mathematica to set up this system of equations and get the error associated with the minimization? Also, the solution needs to be constrained such that C$_n$ is positive.

phosBB = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_BB.CSV"];
phosMR = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_MR.CSV"];
phosPT = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_PT.CSV"];
phosTB = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_TB.CSV"];
phosUI = Import["/home/marco/LatexDocs/analytical/uvvis/PHOS_UI.CSV"];


PreProcess data to get rid of headings:

phosBB = Delete[phosBB, 1];
phosMR = Delete[phosMR, 1];
phosPT = Delete[phosPT, 1];
phosTB = Delete[phosTB, 1];
phosUI = Delete[phosUI, 1];


Calculate $\epsilon$ from the known solutions (BB-TB):

\[Epsilon]PhosBB = phosBB[[All, 2]]/concPhosBB;
\[Epsilon]PhosMR = phosMR[[All, 2]]/concPhosMR;
\[Epsilon]PhosPT = phosPT[[All, 2]]/concPhosPT;
\[Epsilon]PhosTB = phosTB[[All, 2]]/concPhosTB;


Setup data table [$A_{unk} \; \epsilon_1 \; \epsilon_2 \; \epsilon_3 \; \epsilon_4$]

phosDataTable =
Transpose[{phosUI[[All, 1]], \[Epsilon]PhosBB, \[Epsilon]PhosMR, \[Epsilon]PhosPT, \[Epsilon]PhosTB}];

-
Have you attempted something already with e.g. FindFit? Can you show your code? –  Mike Honeychurch Feb 22 '14 at 2:37
It looks like a linear system -- so you can rewrite your equations in matrix form m.x==b and then use LinearSolve. m is your data {A, e1, e2, e3, e4}, b is zero and x={c1, c2, c3, c4}. –  bill s Feb 22 '14 at 3:43
@MikeHoneychurch I'm not sure what mathematica function to use. The idea is to guess C1-C4 to try to get as close to zero as possible. There will be some error associated with the calculation because I don't think there is any exact solution. we just want to get close. –  olliepower Feb 22 '14 at 21:14
@bills I tried using Linear solve where m was {e1,e2,e3,e4}, x = {c1,c2,c3,c4}, b = {A\$_{unk}}. (a rearrangement of the equation i posted above). There was no exact solution. –  olliepower Feb 22 '14 at 21:16
400 points is overkill for visible spectra of four components. If you have access to Harris' Quantitative Chemical Analysis, chapter 18 provides a very thorough explanation and demonstrates the solution using Excel (which turns out to be a MMA one or two liner). –  bobthechemist Feb 23 '14 at 0:21

Import the data and filenames, cull the datasets of information that is not useful:

data = Import /@ FileNames["*.csv"];
filenames = FileNames["*.csv"];
data = data[[All, 2 ;;, {1, 2}]];
ListLinePlot[data, PlotLegends -> filenames]


Looks like PT data was either incorrectly collected or is effectively transparent in the region of interest. The question currently is missing concentration information, so I will assume each standard is 0.1 M just for pedagogical purposes. Convert the standard ordinates to molar absorptivities.

conc = {0.1, 0.1, 0.1, 0.1};
MapIndexed[(data[[#2, All, 2]] /= #1) &, conc];


Now solve the linear system. Here I'm using the definition of the linear lease squares just for fun.

soln = With[{A = data[[1 ;; 4, All, 2]]},
Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data[[5, All, 2]],
{c1, c2, c3, c4}]]
(* {{c1 -> 0.059183, c2 -> 0.00586333, c3 -> 0.0171604, c4 -> 0.0195068}} *)


(* Note, I would have liked to replace {c1, c2, c3, c4} with the actual filenames but they are not valid mathematica symbols *)

The above expression can also (should also?) be entered as LeastSquares[Transpose@data[[1 ;; 4, All, 2]], data[[5, All, 2]]] which yields the same result as a list of coefficients rather than a list of rules.

Create the theoretical curve and compare it to the original

theoretical = Transpose[{data[[1, All, 1]],
{c1, c2, c3, c4}.data[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data[[5]]}]


If you are interested in parameter errors, then use LinearModelFit.

lm = LinearModelFit[{Transpose@data[[1 ;; 4, All, 2]], data[[5, All, 2]]}];
lm["ParameterTable"]


Not surprisingly, the third variable, corresponding to PT, has the largest error.

-
you are right about the PT. The csv file is still saved on the instrument. I couldn't get into the room to retrieve the spectrum. Thank you so much. I'll repost tomorrow with the PT included so that we can get a good answer. –  olliepower Feb 23 '14 at 3:17
I consulted my notebook. The PT solution was clear at at this PH. Could you elaborate on the theoretical and actual fit plot? I'm not sure I understand what is being plotted. Thank you so much. I'll repost tomorrow with all of the data processed in case you are interested. –  olliepower Feb 23 '14 at 3:49
@olliepower the theoretical plot is simulated from the molar absorptivities of the standard compounds and the best fit concentrations of these standards that makes up the final solution. The other line is the experimental plot. –  bobthechemist Feb 23 '14 at 14:30

Here is the final result if you are interested.

 data = Import /@
FileNames["/home/marco/LatexDocs/analytical/uvvis/PHOS_*.CSV"];
data = data[[All, 2 ;;, {1, 2}]];
ListLinePlot[data, PlotLegends -> {"BB", "MR", "PT", "TB", "UI"},
AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16],
Style["Absorbance", FontSize -> 16]},
PlotRange -> {{400, 750}, {0, 0.9}}]
concBB = 1.6*10^(-5);
concMR = 37.1  * 10^(-5);
concPT = 31.4 * 10^(-5);
concTB =  21.4 * 10^(-5);
conc = {concBB, concMR, concPT, concTB};
MapIndexed[(data[[#2, All, 2]] /= #1) &, conc];
soln = With[{A = data[[1 ;; 4, All, 2]]},
Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data[[5, All, 2]], {c1,
c2, c3, c4}]];

theoretical =
Transpose[{data[[1, All,
1]], {c1, c2, c3, c4}.data[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data[[5]]},
AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16],
Style["Absorbance", FontSize -> 16]}]

lm = LinearModelFit[{Transpose@data[[1 ;; 4, All, 2]],
data[[5, All, 2]]}];
error1 = Total[(#^2) & @  lm["FitResiduals"]];

data2 = Import /@
FileNames["/home/marco/LatexDocs/analytical/uvvis/CARB_*.CSV"];
data2 = data2[[All, 2 ;;, {1, 2}]];
ListLinePlot[data2, PlotLegends -> {"BB", "MR", "PT", "TB", "UI"},
AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16],
Style["Absorbance", FontSize -> 16]},
PlotRange -> {{400, 750}, {0, 0.9}}]

MapIndexed[(data2[[#2, All, 2]] /= #1) &, conc];
soln2 = With[{A = data2[[1 ;; 4, All, 2]]},
Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data2[[5, All, 2]], {c1,
c2, c3, c4}]];

theoretical =
Transpose[{data2[[1, All,
1]], {c1, c2, c3, c4}.data2[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data2[[5]]}]
lm2 = LinearModelFit[{Transpose@data2[[1 ;; 4, All, 2]],
data2[[5, All, 2]]}];
error2 = Total[(#^2) & @  lm2["FitResiduals"]];

data3 = Import /@
FileNames["/home/marco/LatexDocs/analytical/uvvis/KHP_*.CSV"];
data3 = data3[[All, 2 ;;, {1, 2}]];
ListLinePlot[data3, PlotLegends -> {"BB", "MR", "PT", "TB", "UI"},
AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16],
Style["Absorbance", FontSize -> 16]},
PlotRange -> {{400, 750}, {0, 0.9}}]

MapIndexed[(data3[[#2, All, 2]] /= #1) &, conc];
soln3 = With[{A = data3[[1 ;; 4, All, 2]]},
Solve[A.Transpose@A.{c1, c2, c3, c4} == A.data3[[5, All, 2]], {c1,
c2, c3, c4}]];

theoretical =
Transpose[{data3[[1, All,
1]], {c1, c2, c3, c4}.data3[[1 ;; 4, All, 2]] /. First@soln}];
ListLinePlot[{theoretical, data3[[5]]},
AxesLabel -> {Style["wavelength (nm)"  , FontSize -> 16],
Style["Absorbance", FontSize -> 16]}]
lm3 = LinearModelFit[{Transpose@data3[[1 ;; 4, All, 2]],
data3[[5, All, 2]]}];
error3 = Total[(#^2) & @  lm3["FitResiduals"]];

dataSummary = {c1, c2, c3, c4} /. {soln, soln2, soln3}
errorSummary = {error1, error2, error3}


Mean Concentration Molar

{{0.000469681, 0.00678397, 0.0169567, 0.0012554}}


Plots: each set of three corresponds to: (spectrum, theoretical, plot of residuals)

-
@bobthechemist I have updated the final result if you are interested. Is there any reason why the theoretical plot looks so good for the first plot and so bad for the rest of them? Its pretty odd that there is an apparent structure to the residuals as well. –  olliepower Feb 23 '14 at 18:31