## dharr

Dr. David Harrington

## 1375 Reputation

15 years, 343 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

## step by step...

Seems there should be a simpler way, but you can do it step by step.

## missing assignment...

sigma[CNT] = 0.1e7; should be sigma[CNT] := 0.1e7;

In general it is helpful for us to solve your problems if you upload your worksheet using the fat green up arrow in the editor.

## more Digits and better guesses...

Since I needed to figure out good guesses and you formulated the ionic strength in molality and the others in moles I got frustrated and cheated by refomulating the whole problem in molalities, which is how I would have set up this problem anyway. These equilibrum problems almost always need high Digits and good guesses. You will see some chemical intuition was also required to solve it.

[Edit: at the end I check it satisfies the OPs eqns]

 > restart;

Given equations rewritten in terms of molalities m_X (=n_X/m). Omit x and just solve directly for the actual concentrations in the solution.

 > eq1:=K_1=(m_Cl*u_Cl)*(m_H*u_H)/(m_HCl):    #K1 for HCl(aq)->H+ + Cl- eq2:=K_2=(m_Na*u_Na)*(m_OH*u_OH)/(m_NaOH): #K2 for NaOH(aq)->Na+ + OH- eq3:=K_w=(m_H*u_H)*(m_OH*u_OH):            #Kw for H2O -> H+ + OH- (neglect water activity effects) eq4:=m_NaCl=m_Na+m_NaOH:                   #mass balance for Na eq5:=m_NaCl=m_Cl+m_HCl:                    #mass balance for Cl eq6:=m_Na+m_H=m_Cl+m_OH:                   #charge balance eq7:=2*ionic=m_H+m_Cl+m_Na+m_OH:           #definition of ionic strength eq8:=u_H=0.4077*ionic^2-0.3152*ionic+0.9213:  #activity coeff of H+ eq9:=u_Na=0.0615*ionic^2-0.2196*ionic+0.8627: #activity coeff of Na+ eq10:=u_OH=0.1948*ionic^2-0.1803*ionic+0.8887:#activity coeff of HO- eq11:=m=r*V:                                  #estimate of mass of water from density and volume eq12:=u_Cl=(1.417625986641341e-01)*exp(-ionic/2.199955601666953e-02)+ #activity coeff of Cl-       2.369460669647978e-01*exp(-ionic/3.756472377688394e-01)+5.859738096037875e-01:

Given parameters

 > params:={K_1=1.*10^8,K_2=10^(-0.2),K_w=10^(-13.995),n_NaCl=1.2*10^(-4),V=5.*10^(-8),r=997.0,x=0}:

m is directly available from these parameters and eq11.

 > {m=eval(rhs(eq11),params)};params:=params union %: And with m we can get the initial molality of NaCl

 > {m_NaCl=eval(n_NaCl/m,params)};params:=params union %: What molalities are we actually interested in? Na+ and Cl- will be close to the above value, OH- and H+ will be smaller

 > mmax:=eval(m_NaCl,params); So in this concentration range, what are the activity coefficients (u) of H+ (red), Na+ (blue) , OH- (green), Cl- (magenta)

 > plot(map(rhs,[eq8,eq9,eq10,eq12]),ionic=0..3,color=[red,blue,green,magenta]); Leave eq 11 out since we know m. We will specify generous, but not too generous, ranges for fsolve

For the activities, the above plots suggest 0..4, but since thay are monotonic and we know we are above the minimum it may be safer to specify a restricted range for "ionic"

We choose 2*mmax for Na+ and Cl-  molalities (0..nmax should work but we don't want to be too close to the end of the range) , and 0.5*mmax...2*mmax for ionic (could choose this also for Na+ and Cl-)

For the OH- , H+, HCl, and maybe NaOH concentrations we can presumably go much smaller (msmall)

 > nmax:=2*eval(mmax*m,params); > msmall:=1e-2; Since the equilibrium constants and other parameters vary over many orders of magnitude, this should be solved with a high setting of Digits

 > Digits:=25; > eqns:=eval({eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq12},params):

First attempt gives no solution

 > res:=fsolve(eqns,{m_HCl=0..msmall,m_Cl=0..mmax,m_H=0..msmall,m_Na=0..mmax,m_OH=0..msmall,m_NaOH=0..mmax,ionic=0.5*mmax..2*mmax,u_H=0..4,u_Na=0..4,u_OH=0..4,u_Cl=0..4}); Solve the simpler problem by assuming that there is no HCl(aq) in solution - so omit eq 1 and set m_HCl=0 in eq 5

 > eqns1:=eval({eq2,eq3,eq4,eval(eq5,m_HCl=0),eq6,eq7,eq8,eq9,eq10,eq12},params):
 > res1:=fsolve(eqns1,{m_Cl=0..mmax,m_H=0..msmall,m_Na=0..mmax,m_OH=0..msmall,m_NaOH=0..mmax,ionic=0.5*mmax..2*mmax,u_H=0..4,u_Na=0..4,u_OH=0..4,u_Cl=0..4}); Values look reasonable. Can we use this as a starting guess for the problem with HCl - guess m_HCl=0

 > guess:={m_HCl=0} union res1:
 > res2:=fsolve(eqns,guess); Looks good - HCl is very small and the change in m_Cl shows the requirement to use many digits. (won't work with Digits:=10)

Convert to the original mol units of the OP, and check it solves the OPs equations

 > molans:=eval({n_HCl=m_HCl*m,n_H=m_H*m,n_Cl=m_Cl*m,n_Na=m_Na*m,n_NaOH=m_NaOH*m,n_OH=m_OH*m},params union res2) union res2; > eeq1:=K_1=(((n_Cl-x)*u_Cl)*(n_H*u_H))/(n_HCl*m): eeq2:=K_2=(((n_Na-x)*u_Na)*(n_OH*u_OH))/(n_NaOH*m): eeq3:=K_w=(n_H*u_H/m)*(n_OH*u_OH/m): eeq4:=(n_NaCl-x)=(n_Na-x)+n_NaOH: eeq5:=(n_NaCl-x)=(n_Cl-x)+n_HCl: eeq6:=(n_Na-x)+n_H=(n_Cl-x)+n_OH: eeq7:=2*ionic=(n_H/m)+((n_Cl-x)/m)+((n_Na-x)/m)+(n_OH/m): eeq8:=u_H=0.4077*ionic^2-0.3152*ionic+0.9213: eeq9:=u_Na=0.0615*ionic^2-0.2196*ionic+0.8627: eeq10:=u_OH=0.1948*ionic^2-0.1803*ionic+0.8887: eeq11:=m=r*V: eeq12:=u_Cl=(1.417625986641341e-01)*exp(-ionic/2.199955601666953e-02)+2.369460669647978e-01*exp(-ionic/3.756472377688394e-01)+5.859738096037875e-01:

We find that it does.

 > OPans:=eval([eeq1,eeq2,eeq3,eeq4,eeq5,eeq6,eeq7,eeq8,eeq9,eeq10,eeq11,eeq12],molans union params); >

## reduction by n!...

There is no difference in the singularity if the columns are permuted, so choosing the entries of the first row without considering different orders reduces the number of ways by n! for an n x n matrix.

 > restart:
 > alias(det = LinearAlgebra[Determinant]); with(combinat): Only test up to column permutations. Find all combinations to fill the first row, and then permutations of the other entries. Reduces the number of ways by a factor of n! for an nxn Matrix

 > st:=time(): S    := 0: n    := 3; N    := n^2: COMB:=choose(N,n): nums:=[\$1..N]: for comb in COMB do   plist:=remove(has,nums,comb);   PERM:=permute(plist);   for perm in PERM do     M:=Matrix(n,n,ListTools:-Flatten([comb,perm]));     if det(M)=0 then S := S+1; end if:   end do: end do: S*n!; time()-st;   Brute Force

 > st:=time(): S    := 0: n    := 3; N    := n^2: PERM := permute(N): for perm in PERM do   M := Matrix(n, n, perm):   if det(M)=0 then S := S+1; end if: end do: S; time()-st;   >

## submatrix...

T is already a Matrix, and stripping off the first two rows and columns can be done with T[3..,3..]                 ## normalization=full...

You need normalization=full, and then, for some reason, another factor of 2 (Every time I do a DFT in a different piece of software, I always need to do a sample to figure this out).

FFT.mw

## Should work...

I couldn't see why yours didn't work, but here's an example that does. The button code is

```use DocumentTools in
nrows:=nrows+1;
A:=Matrix(nrows,2,A);
A[nrows,1]:=0;
A[nrows,2]:=1;
Do(%DataTable0(VisibleRows)=nrows);
end use;
```

## force floating point...

If your change ln(2) to ln(2.) then that will force the answer into a single numerical floating point value.

## Document properties...

In the file menu under doccument properties, you need an attibute Active, which has the value true for appearance in the help browser as a worksheet, and false for appearance as a help page.

## with(plots)...

You need with(plots): since implicitplot, pointplot and display are all in the plots package. Using with(plots,implicitplot) loads only implicitplot and not the others.

## 0 or 1...

Originally _B1, renamed _B1~:
is assumed to be: OrProp(0,1)

so it means either zero or 1

## D...

 > restart;
 > with(plots): with(VectorCalculus): rf:=t-> < cos(t), sin(t), t >; > rCurve := spacecurve( rf(t), t = -1 .. 1 ); > drf:=D(rf):
 > drCurve := spacecurve( drf(t), t = -1 .. 1 ); >

## addcoords only works one way...

to use addcoordinates you need to know x, y and z in terms of the new coordinate system, it doesn't do inverse ones. pdechangecoords is deprecated, but dchange can do what you want, though the form of the inverse transformation is a mess:

 > restart;
 > with(PDEtools):
 > toroidaltrans:={x=a*sinh(v)*cos(w)/(cosh(v)-cos(u)),y=a*sinh(v)*sin(w)/(cosh(v)-cos(u)),z=a*sin(u)/(cosh(v)-cos(u))}; Solving tor the toroidal oordinates in terms of the cartesian coordinates gives a mess. Adding assumptions and using simplify with options might help here.

 > invtoroidaltrans:=solve(toroidaltrans,{u,v,w}); > fn:=diff(f(x,y,z),z); To toroidal

 > fntoroidal:=dchange(toroidaltrans,fn,[u,v,w],params={a}); And back again

 > fnback:=dchange(invtoroidaltrans,fntoroidal,[x,y,z],itr=toroidaltrans,params={a},simplify); >

## randomly multiply by matrices...

Perhaps something like this...

 > restart;
 > R:=Matrix(3,3,[1,0,0,0,cos(delta),sin(delta),0,-sin(delta),cos(delta)]); R:=Matrix(3,3,[cos(delta),0,-sin(delta),0,1,0,sin(delta),0,cos(delta)]);  Random number generators to select which matrix and which delta value

 > r01:=rand(0..1):r2Pi:=rand(0..2.*Pi):

Successively multiply by a randomly selected Matrix with a randomly selected delta

 > N:=1000; S:=table():S:=Vector([0,1,0]); #initial vector for i from 2 to N do S[i]:=eval(R[r01()],delta=r2Pi()).S[i-1] end do: S:=convert(S,list):  > plots:-pointplot3d(S,symbolsize=10,color=black); >

Edit: altered to make initial vector a unit vector.

I came to the same conclusion as @vv, that there is not a unique solution. Here is a solution to the simpler problem in which P6 is absent and length L25 is known.

 > restart;

Mapleprimes problem proposed by @Fancypants.

Set up some viable points just to figure out some viable lengths. Make a drawing of the problem.
Problem: the lengths of the lines are given, as are P0 and P1. Find the locations of the other points.

 > with(geometry): pxy:={p3x=-0.5,p3y=3,p5x=3,p5y=4,p6x=3.5,p6y=3}: point(P0,0,0):point(P1,1,0):point(P3,eval([p3x,p3y],pxy)):point(P5,eval([p5x,p5y],pxy)):point(P6,eval([p6x,p6y],pxy)): OnSegment(P2,P0,P3,1.2):OnSegment(P4,P3,P5,0.7): segment(P0P1,P0,P1):segment(P0P3,P0,P3):segment(P3P5,P3,P5):segment(P5P6,P5,P6):segment(P1P4,P1,P4):segment(P2P6,P2,P6): plots:-display(draw([P0,P1,P2,P3,P4,P5,P6],printtext=true,axes=none),                draw([P0P1,P0P3,P3P5,P5P6,P1P4,P2P6],axes=none)); lengths=[L01=distance(P0,P1),L02=distance(P0,P2),L23=distance(P2,P3),L34=distance(P3,P4),L45=distance(P4,P5),L56=distance(P5,P6),L26=distance(P2,P6),L14=distance(P1,P4),L25=distance(P2,P5)];  Now solve it only knowing these lengths.

 > restart;
 > with(geometry):
 > assume(L01>0,L02>0,L23>0,L34>0,L45>0,L14>0,L56>0,L26>0,L25>0,p3y>0,p5y>0,p6y>0);

Given lengths

 > lengths:={L01 = 1, L02 = 1.658935235, L23 = 1.382446030, L34 = 1.498846154, L45 = 2.141208791, L56 = 1.118033989, L26 = 4.011605067, L14 = 3.412271768, L25 = 4.037018784}: #assign(lengths);

P0 and P1 are known, set P0 at the origin and P1 at (L01,0)

 > point(P0,0,0);point(P1,L01,0);  The points we need to find

 > point(P3,p3x,p3y);point(P5,p5x,p5y);point(P6,p6x,p6y);   Now define some segments we need to draw for the two quadrilaterals (the segments are just for drawing, not for calculation).

 > segment(P0P1,P0,P1);segment(P0P3,P0,P3);segment(P3P5,P3,P5);segment(P5P6,P5,P6);    Point P2 is on line P0P3 at ratio k = L02/L23 and similarly for point P4

 > OnSegment(P2,P0,P3,L02/L23); OnSegment(P4,P3,P5,L34/L45);  Now we can define the other two segments

 > segment(P1P4,P1,P4);segment(P2P6,P2,P6);  Setup the equations to solve for the locations of all the points. P0 and P1 are known immediately in terms of L01, and P2 and P4 are given in terms of the others, so  the 6 unknowns are the coordinates of P3, P5 and P6.

Specifying all the lengths gives only 5 independent equations - need another

 > eq02:=L02^2=simplify(distance(P0,P2)^2); eq23:=L23^2=simplify(distance(P2,P3)^2); #equiv to eq02 eq34:=L34^2=simplify(distance(P3,P4)^2); eq45:=L45^2=simplify(distance(P4,P5)^2); #equiv to eq34 eq14:=L14^2=simplify(distance(P1,P4)^2); eq56:=L56^2=simplify(distance(P5,P6)^2); eq26:=L26^2=simplify(distance(P2,P6)^2); eq25:=L25^2=simplify(distance(P2,P5)^2); #not supposed to know this length        First solve the slightly simpler problem, where we know the length L25 but don't have point P6 in the problem. Now we have 4 equations and 4 unkowns

 > eqs:={eq02,eq34,eq14,eq25}:nops(eqs); indets(eqs,name);  Solve takes forever, so try fsolve, which returns the correct solution

 > fsolve(eval(eqs,lengths),{p3x=-2..0,p3y=0..10,p5x=0..10,p5y=0..10}); Is there a unique solution with the information given, i.e. is the system rigid?

Adding the point P6 and replacing P2P5 by two other sides of the triangle P2P5P6 seems to allow for some flexibility?

 1 2 3 4 5 6 7 Last Page 1 of 26