Advanced Search

dharr

Dr. David Harrington

1375 Reputation

13 Badges

15 years, 348 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
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.

MaplePrimes Activity


These are answers submitted by dharr

step by step...

dharr 1375
July 29 2020
1 0

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

Download sums.mw

missing assignment...

dharr 1375
July 28 2020
1 2

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...

dharr 1375
July 26 2020
2 0

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 %:

{m = 0.4985000000e-4}

And with m we can get the initial molality of NaCl

> 

{m_NaCl=eval(n_NaCl/m,params)};params:=params union %:

{m_NaCl = 2.407221665}

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);

2.407221665

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);

0.2400000000e-3

> 

msmall:=1e-2;

0.1e-1

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;

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});

fsolve({0.1011579454e-13 = m_H*u_H*m_OH*u_OH, .6309573445 = m_Na*u_Na*m_OH*u_OH/m_NaOH, 2.407221665 = m_Cl+m_HCl, 2.407221665 = m_Na+m_NaOH, 100000000. = m_Cl*u_Cl*m_H*u_H/m_HCl, u_Cl = .1417625986641341*exp(-45.45546279*ionic)+.2369460669647978*exp(-2.662072017*ionic)+.5859738096037875, u_H = .4077*ionic^2-.3152*ionic+.9213, u_Na = 0.615e-1*ionic^2-.2196*ionic+.8627, u_OH = .1948*ionic^2-.1803*ionic+.8887, 2*ionic = m_H+m_Cl+m_Na+m_OH, m_Na+m_H = m_Cl+m_OH}, {ionic, m_Cl, m_H, m_HCl, m_Na, m_NaOH, m_OH, u_Cl, u_H, u_Na, u_OH}, {ionic = 1.2036108325 .. 4.814443330, m_Cl = 0 .. 2.407221665, m_H = 0 .. 0.1e-1, m_HCl = 0 .. 0.1e-1, m_Na = 0 .. 2.407221665, m_NaOH = 0 .. 2.407221665, m_OH = 0 .. 0.1e-1, u_Cl = 0 .. 4, u_H = 0 .. 4, u_Na = 0 .. 4, u_OH = 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});

{ionic = 2.407221687118782060747778, m_Cl = 2.407221665000000000000000, m_H = 0.1143810110208346910486800e-6, m_Na = 2.407221572737771039913087, m_NaOH = 0.9226222896008691260626429e-7, m_OH = 0.2211878206074777844241567e-7, u_Cl = .5863642949085324404566185, u_H = 2.525049539726357549337178, u_Na = .6904491669412176811311638, u_OH = 1.583488655494620712846916}

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);

{ionic = 2.407221687118778378231341, m_Cl = 2.407221664999995923312451, m_H = 0.1143810089824910411470732e-6, m_HCl = 0.4076687549288832532913356e-14, m_Na = 2.407221572737769395740300, m_NaOH = 0.9226223060425970007002799e-7, m_OH = 0.2211878245491889036587774e-7, u_Cl = .5863642949085324442845943, u_H = 2.525049539726351481844259, u_Na = .6904491669412173994618614, u_OH = 1.583488655494617923143445}

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;

{ionic = 2.407221687118778378231341, m_Cl = 2.407221664999995923312451, m_H = 0.1143810089824910411470732e-6, m_HCl = 0.4076687549288832532913356e-14, m_Na = 2.407221572737769395740300, m_NaOH = 0.9226223060425970007002799e-7, m_OH = 0.2211878245491889036587774e-7, n_Cl = 0.1200000000002497967771257e-3, n_H = 0.5701893297777178401181599e-11, n_HCl = 0.2032228743320483017657308e-18, n_Na = 0.1199999954009778043776540e-3, n_NaOH = 0.4599272195622346048490895e-11, n_OH = 0.1102621305377706684739005e-11, u_Cl = .5863642949085324442845943, u_H = 2.525049539726351481844259, u_Na = .6904491669412173994618614, u_OH = 1.583488655494617923143445}

> 

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);

[100000000. = 100000000.0000000000000000, .6309573445 = .6309573445000000000000000, 0.1011579454e-13 = 0.1011579454000000000000000e-13, 0.1200000000e-3 = 0.1200000000002500000000000e-3, 0.1200000000e-3 = 0.1200000000002500000000000e-3, 0.1200000011028711021548324e-3 = 0.1200000011028711021548324e-3, 4.814443374237556756462682 = 4.814443374237556756462683, 2.525049539726351481844259 = 2.525049539726351481844259, .6904491669412173994618614 = .6904491669412173994618614, 1.583488655494617923143445 = 1.583488655494617923143445, 0.4985000000e-4 = 0.49850000000000e-4, .5863642949085324442845943 = .5863642949083517098635012]

> 

 

Download Chem.mw

reduction by n!...

dharr 1375
July 22 2020
1 6

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):

det

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;

3

2736

4.672

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;

3

2736

32.110

> 

 

 

Download HowManyMatricesAreSingular.mw

submatrix...

dharr 1375
July 21 2020
1 8

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

with(DifferentialGeometry); with(LieAlgebras)

NULL

M := [Matrix([[1, 0], [0, 0]]), Matrix([[0, 0], [0, 1]]), Matrix([[0, 1], [0, 0]])]

M := [Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = 0, (1, 2) = 1, (2, 1) = 0, (2, 2) = 0})]

L := LieAlgebraData(M, Alg1)

_DG([["LieAlgebra", Alg1, [3, table( [ ] )]], [[[1, 3, 3], 1], [[2, 3, 3], -1]]])

DGsetup(L)

`Lie algebra: Alg1`

T := MultiplicationTable("LieTable")

T := Matrix(5, 5, {(1, 1) = ``, (1, 2) = `| `, (1, 3) = e1, (1, 4) = e2, (1, 5) = e3, (2, 1) = ``, (2, 2) = `----`, (2, 3) = `----`, (2, 4) = `----`, (2, 5) = `----`, (3, 1) = e1, (3, 2) = `| `, (3, 3) = 0, (3, 4) = 0, (3, 5) = Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (4, 1) = e2, (4, 2) = `| `, (4, 3) = 0, (4, 4) = 0, (4, 5) = Typesetting:-mcomplete(Typesetting:-mrow(Typesetting:-mo("&uminus0;"), Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DGproduct"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mi("e3"))))), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("&uminus0;1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (5, 1) = e3, (5, 2) = `| `, (5, 3) = Typesetting:-mcomplete(Typesetting:-mrow(Typesetting:-mo("&uminus0;"), Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DGproduct"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mi("e3"))))), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("&uminus0;1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (5, 4) = Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (5, 5) = 0})

whattype(T)

Matrix

A := T[3 .. (), 3 .. ()]

A := Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (2, 1) = 0, (2, 2) = 0, (2, 3) = Typesetting:-mcomplete(Typesetting:-mrow(Typesetting:-mo("&uminus0;"), Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DGproduct"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mi("e3"))))), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("&uminus0;1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (3, 1) = Typesetting:-mcomplete(Typesetting:-mrow(Typesetting:-mo("&uminus0;"), Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DGproduct"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mi("e3"))))), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("&uminus0;1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (3, 2) = Typesetting:-mcomplete(Typesetting:-mi("e3"), Typesetting:-mrow(Typesetting:-mi("_DG"), Typesetting:-mo("⁡"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-ms("vector"), Typesetting:-mi("Alg1"), Typesetting:-mfenced(open = "[", close = "]"), open = "[", close = "]"), Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mn("3"), open = "[", close = "]"), Typesetting:-mn("1"), open = "[", close = "]"), open = "[", close = "]"), open = "[", close = "]")))), (3, 3) = 0})

LinearAlgebra:-Rank(A)

2

NULL

 

Download CommutatorExample.mw

normalization=full...

dharr 1375
July 19 2020
2 3

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...

dharr 1375
July 16 2020
0 2

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;

Download Morerows.mw

force floating point...

dharr 1375
July 16 2020
0 0

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

Document properties...

dharr 1375
July 15 2020
2 1

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)...

dharr 1375
July 11 2020
1 1

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...

dharr 1375
July 11 2020
2 0

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

so it means either zero or 1

D...

dharr 1375
July 10 2020
2 0


 

> 

restart;

> 

with(plots): with(VectorCalculus):
rf:=t-> < cos(t), sin(t), t >;

proc (t) options operator, arrow; VectorCalculus:-`<,>`(cos(t), sin(t), t) end proc

> 

rCurve := spacecurve( rf(t), t = -1 .. 1 );

> 

drf:=D(rf):
> 

drCurve := spacecurve( drf(t), t = -1 .. 1 );

> 

 

 

Download D.mw

addcoords only works one way...

dharr 1375
July 09 2020
3 0

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))};

{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});

{u = arctan(-2*a*z/((2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2+2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)*RootOf((-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2+a^2+x^2+y^2+z^2)*_Z^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)), (a^2-x^2-y^2-z^2)/((2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2+2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)*RootOf((-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2+a^2+x^2+y^2+z^2)*_Z^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2))), v = ln(RootOf((-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2+a^2+x^2+y^2+z^2)*_Z^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*x^2-2*RootOf((x^2+y^2)*_Z^2-1)*a*y^2-a^2-x^2-y^2-z^2)), w = arctan(y*RootOf((x^2+y^2)*_Z^2-1), RootOf((x^2+y^2)*_Z^2-1)*x)}

> 

fn:=diff(f(x,y,z),z);

diff(f(x, y, z), z)

To toroidal

> 

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

(-(diff(f(u, v, w), u))+(diff(f(u, v, w), u))*cosh(v)*cos(u)-(diff(f(u, v, w), v))*sinh(v)*sin(u))/a

And back again

> 

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

diff(f(x, y, z), z)

> 

 


 

Download toroidal.mw

randomly multiply by matrices...

dharr 1375
July 07 2020
2 0

Perhaps something like this...

> 

restart;
> 

R[0]:=Matrix(3,3,[1,0,0,0,cos(delta),sin(delta),0,-sin(delta),cos(delta)]);
R[1]:=Matrix(3,3,[cos(delta),0,-sin(delta),0,1,0,sin(delta),0,cos(delta)]);

R[0] := Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = cos(delta), (2, 3) = sin(delta), (3, 1) = 0, (3, 2) = -sin(delta), (3, 3) = cos(delta)})

R[1] := Matrix(3, 3, {(1, 1) = cos(delta), (1, 2) = 0, (1, 3) = -sin(delta), (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = sin(delta), (3, 2) = 0, (3, 3) = 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[1]:=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):

N := 1000

S[1] := Vector(3, {(1) = 0, (2) = 1, (3) = 0})

> 

plots:-pointplot3d(S,symbolsize=10,color=black);

> 

 


Edit: altered to make initial vector a unit vector.

Download Stokes.mw

partial answer...

dharr 1375
July 07 2020
0 1

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)];

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]

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);

P0

P1

The points we need to find

> 

point(P3,p3x,p3y);point(P5,p5x,p5y);point(P6,p6x,p6y);

P3

P5

P6

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);

P0P1

P0P3

P3P5

P5P6

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);

P2

P4

Now we can define the other two segments

> 

segment(P1P4,P1,P4);segment(P2P6,P2,P6);

P1P4

P2P6

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

L02^2 = L02^2*(p3x^2+p3y^2)/(L23+L02)^2

L23^2 = L23^2*(p3x^2+p3y^2)/(L23+L02)^2

L34^2 = (p3y^2-2*p3y*p5y+p5y^2+(p3x-p5x)^2)*L34^2/(L45+L34)^2

L45^2 = L45^2*(p3y^2-2*p3y*p5y+p5y^2+(p3x-p5x)^2)/(L45+L34)^2

L14^2 = ((L01^2-2*L01*p5x+p5x^2+p5y^2)*L34^2+2*L45*(L01^2+(-p3x-p5x)*L01+p3x*p5x+p3y*p5y)*L34+L45^2*(L01^2-2*L01*p3x+p3x^2+p3y^2))/(L45+L34)^2

L56^2 = p5y^2-2*p5y*p6y+p6y^2+(p5x-p6x)^2

L26^2 = (L02^2*(p3x^2-2*p3x*p6x+p3y^2-2*p3y*p6y+p6x^2+p6y^2)-2*(-p6y^2+p3y*p6y+p6x*(p3x-p6x))*L23*L02+L23^2*(p6x^2+p6y^2))/(L23+L02)^2

L25^2 = (L02^2*(p3x^2-2*p3x*p5x+p3y^2-2*p3y*p5y+p5x^2+p5y^2)-2*(-p5y^2+p3y*p5y+p5x*(p3x-p5x))*L23*L02+L23^2*(p5x^2+p5y^2))/(L23+L02)^2

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);

4

{L01, L02, L14, L23, L25, L34, L45, p3x, p3y, p5x, p5y}

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});

{p3x = -.4999999995, p3y = 3.000000000, p5x = 2.999999999, p5y = 4.000000002}

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?

 

Download geometry2.mw

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