Errors in Addition and Subtraction
Hi All,
I'm attempting to numerically model some electrochemical equilibria but am getting bizarre and inconsistent output from my code. The basic idea is that I start with j species (in my case, 4 species, so j ranges from 0 to 3) with populations {1,0,0,0} in the first row of a two-dimensional wave I've named Species with rows 0 to i. I use an equation to calculate what the change in population of each species should be in the current row based on the contents of the previous row. For testing purposes, I've written the delta between rows to a separate wave that I've named Delta. I then add the delta from the previous row to the current one and subtract the delta from the previous row and column to get the population of the current row of the previous column.
The problem is - the subtraction occurs but the addition does not. Moreover, when the code moves from Species[1][1] to Species[1][2], it subtracts the Delta[1][1] value from Species[1][1] and not the Delta[1][2] value (which is 0). The net result is that Species[*][0] decays as expected, but all of the other columns just have values of 0 for every entry. I integrated the Delta[*][1] values and they give the desired rising sigmoidal curve I'd expected for Species[*][1], but Species[*][1] itself is just 0.
If you'd like to try to reproduce this, the code relies on a wave "Ered" which is {-0.38,-0.61,-0.72} that I made manually.
Function MakePopWaves(ERed) //Automatically populate waves for population based on a set of reduction potentials Wave ERed Variable NumERed = numpnts(ERed) Variable NumSpecies = NumERed + 1 Variable i Make/D/O/N=(1001,NumSpecies-1) KVals Make/D/O/N=(1001,NumSpecies) Species SetScale/p x 0,-0.001,"",KVals SetScale/p x 0,-0.001,"",Species String PathToERed = GetWavesDataFolder(ERed,2) String PathToKVals = GetWavesDataFolder(KVals,2) String PathToSpecies = GetWavesDataFolder(Species,2) PopModel(PathToERed,PathToKVals,PathToSpecies) End Function PopModel(PathToERed,PathToKVals,PathToSpecies) String PathToERed,PathToKVals,PathToSpecies WAVE ERed = $PathToERed WAVE KVals = $PathToKVals WAVE Species = $PathToSpecies Species = 0 Species[0][0] = 1 Duplicate/D/O Species,Delta Delta = 0 Variable NumERed = numpnts(ERed) Variable WaveLen = DimSize(Species,0) Variable NumSpecies = NumERed + 1 Variable i,j Variable PotDelta = DimDelta(Species,0) Variable PotOffset = DimOffset(Species,0) Make/D/O/FREE/N=(WaveLen) Potentials SetScale/p x PotOffset,PotDelta,"",Potentials Potentials = x Variable R = 8.314 //J/mol K Variable T = 298 // K Variable z = 1 // electron transferred Variable F = 96485 // C/mol Display; Variable/D Nernst //The equation is getting unwieldy, calc a portion of the Nernst equation separately to make it easier For(i=1;i<WaveLen;i+=1) For(j=1;j<NumERed;j+=1) Nernst = exp((ERed[j]-Potentials[i])/(R*T/(z*F))) Delta[i][j] = (Nernst*Species[i-1][j-1]-Species[i-1][j])/(1+Nernst) Species[i][j] = Species[i-1][j] + Delta[i][j] Species[i][j-1] = Species[i-1][j-1] - Delta[i][j] EndFor EndFor End
Thanks for any insights in advance!
Will
A few thoughts to trim your code.
* Move the constants R(=8.31446261815324!!!), TK(=298.15!!!!) to Static Constants at the top of the code block.
* You need not pass the path to wave, just the wave. A better approach as needed is also to make use of SetDataFolder to move into and back out of a data folder, keeping the clutter within the location for the calculations.
* Vectorize to create a matrix for Nernst in advance
variable RTKzF = R*TK/(z*F)
MNerst[][] = exp((ERed[q] - Potentials[p])/RTKzF)
* Do you really need to create a matrix for Delta[i][j] as opposed to just a variable Delta? Following this, you would write: Species[i][j] = Species[i-1][j] + Delta and Species[[i][j-1] - Delta
* Should you do the subtraction before you do the addition? ????
Finally, do you know how to translate explicit for-loops to vectorized implicit approaches?
For example, Species[1,inf][] = Species[p-1][q] + Delta
April 1, 2025 at 12:23 pm - Permalink
What strikes me is that Delta[i][j] is zero in every other iteration of the inner loop. Is this supposed to be like this?
If you change the start values Species[0][1] and Species[0][2] to non zero values (i.e. 1e-10) something happens, although I don't know if the output make sense.
April 2, 2025 at 03:40 am - Permalink
I think the main problem in your code is here:
Species[i][j] will be overwritten in the next j loop. I don't know your model but I guess the total number od species at each i should be constant (mass conservation law or similar), is that correct ?. This is not the case as can be easily verified adding the following instructions in the loop:
(Advice: always check that your code respects conservation laws!) To respect the conservation law I suggest something like :
If you think a bit you will see that this is not quite the same algorithm. Now I don't know much about your model so I can't tell if it is correct but at least the total number of species is conserved now.
April 2, 2025 at 10:14 am - Permalink
Thank you all for the the helpful comments! The code is cleaner now, I learned about static constants and vectorized implicit loops, and Pawel was right - the code was overwriting the results I was trying to generate. Using += and -= has fixed the issue, and I now get the expected population curves.
Best,
Will
April 3, 2025 at 07:39 am - Permalink