I was seeing some oddities when I solved, so I tried solving two different ways:
R=8.3144;
T=298.15;
F=96487;
eps0= 8.854*10^(-12);
epsr=80;
A=1.4*10^(-20); %hamaker constant
a=200*10^(-9); %radius
cpzc=10^(3)*10^(-4); %Cpdi pzc
cpdi=10^(3)*10^(-2); %Cpdi
syms C h
psip=((R*T)/F)*log(cpdi/cpzc);
k=F*sqrt((2*C)/(eps0*epsr*R*T));
Vt= -((A*a)/(12*h))+((64*pi*a*C*R*T)/(k^2))*(tanh((F*psip)/R*T))^2*exp(-C*h);
Vt1= diff(Vt,h,1);
kvals = -1:1;
nk = length(kvals);
sol_h1 = zeros(1,nk, 'sym');
sol_C1 = zeros(1,nk, 'sym');
sol_h1_raw = solve(Vt, h, 'returnconditions', true); %will be in terms of lambertw
k = sol_h1_raw.parameters; %lambertw parameter
for kidx = 1 : nk
K = kvals(kidx);
h_k = subs(sol_h1_raw.h, k, K);
solc = solve( subs(Vt1, h, h_k), C); %do not need return conditions
if isempty(solc)
fprintf('no C for k = %d\n', K);
sol_C1(kidx) = nan;
sol_h1(kidx) = nan;
else
if ~isscalar(solc)
fprintf('%d C solutions for k = %d\n', length(solc), K);
end
sol_C1(kidx) = solc(1);
sol_h1(kidx) = subs(h_k, C, sol_C1(kidx));
end
end
vpa(sol_C1)
vpa(sol_h1)
subplot(2,2,1); plot(kvals, sol_C1, '*-'); title('C vs k, h first');
subplot(2,2,2); plot(kvals, sol_h1, '*-'); title('h vs k, h first');
sol_h2 = zeros(1,nk, 'sym');
sol_C2 = zeros(1,nk, 'sym');
sol_C2_raw = solve(Vt, C, 'returnconditions', true); %will be in terms of log pi*k*2i
k = sol_C2_raw.parameters; %pi*k parameter
for kidx = 1 : nk
K = kvals(kidx);
C_k = subs(sol_C2_raw.C, k, K);
solh = solve( subs(Vt1, C, C_k), h); %do not need return conditions
if isempty(solh)
fprintf('no h for k = %d\n', K);
sol_C2(kidx) = nan;
sol_h2(kidx) = nan;
else
if ~isscalar(solh)
fprintf('%d h solutions for k = %d\n', length(solh), K);
end
sol_h2(kidx) = solh(1);
sol_C2(kidx) = subs(C_k, h, sol_h2(kidx));
end
end
vpa(sol_C2)
vpa(sol_h2)
subplot(2,2,3); plot(kvals, sol_C2, '*-'); title('C vs k, C first');
subplot(2,2,4); plot(kvals, sol_h2, '*-'); title('h vs k, C first');
The non-nan solutions generated here give exact zeros when substituted in to the equations. They are true solutions.
Solving for h first gives a solution involving lambertw(k, expression) . lambertw() is a special function that deals with solutions to x*exp(x) == constant . https://en.wikipedia.org/wiki/Lambert_W_function . When you deal with real-values only, then for some constants there are no real solutions, and for some constants there is one real solution, and for some constants there are two real solutions. The real solutions can only occur for a parameter usually named k (the "branch") being either -1 or 0.
In most cases, when there are two real solutions, the values are quite different. But in this particular case, if you look very closely at the sol_C1 and sol_h1 you will see that the mathematics works out so that the two solutions are identical -- so although solving for h first looks like it gives two solutions, really the two solutions are the same!
Solving for C first gives a solution involving log() of an expression involving -pi*k*2i . In this particular set of expressions, then only k that leads to a real value is k = 0, so you only get one solution.
If you compare sol_C1(2) to sol_C2(2) you will see that the solutions are identical -- so solving first for h gives the same result as solving first for C. Which is what you would hope for, really.
So, why was I having problems? Because of this:
trysol = vpasolve([Vt,Vt1], [h,C], [1e5; 1e-4]); vpa([trysol.h, trysol.C])
At first look, those appear to be an additional solution "close" to the other solutions. Like 0.00000000006747114290204450599689804038282 vs 0.0000000000000019499416869170421695503665175389 are not so different, right? Maybe my solve() approach missed some solutions, perhaps?
However, if you look more carefully, in the solve() approaches, h is small and C is large, but in this vpasolve() solution, h is large and C is small ! Very different solutions !!
What went wrong? How did vpasolve() find a solution that going forward and backwards using solve() claims is not there? Is solve() missing something?
Let us cross-check on the numeric solution:
vpa(subs([Vt, Vt1], trysol))
Practically zero, looks like there is just some round-off! This solution looks pretty good!
... Except it isn't a real solution. The functions do not cross zero there. What is happening is that vpasolve() works to a numeric tolerance (thinking that it is just dealing with numeric round-off), and as a result, vpasolve() will say a solution exists when the values in the expression get "close enough" to zero. vpasolve() does not prove that the expression mathematically crosses zero.
It is like considering solutions to x^2 + delta for positive delta. If delta is small enough, then any numeric approach will say "oh, look, the function value is very small here, we must just be dealing with some round-off error, there is a solution here (e.g., at 0) to within our tolerance". And yet mathematically there is no real-valued solution to x^2 + delta for positive delta.
vpasolve() is finding a solution within mathematical tolerance given its approach. It just so happens that this is a situation where there is no actual solution in that area.
But
trysol2 = vpasolve([Vt,Vt1], [h,C], [1e-10; 1e10]); vpa([trysol.h, trysol.C]); vpa([trysol2.h, trysol2.C]), vpa(subs([Vt, Vt1], trysol2))
Those residue outputs are not so different from the trysol outputs... but it just so happens that these locations are pretty close to the locations where there are actual mathematical zeros.
How do you know? Unfortunately, without a bunch of plotting or mathematical solutions, you do not know. :( :(
