Stiff ODE: Modified Poisson-Boltzmann
- To: mathgroup at smc.vnet.net
- Subject: [mg118686] Stiff ODE: Modified Poisson-Boltzmann
- From: Brian Giera <bgiera at gmail.com>
- Date: Sat, 7 May 2011 07:32:48 -0400 (EDT)
Physical picture:
I am attempting to solve a set of equations that will give me the voltage
and concentration profiles of hard-sphere (Carnahan-Starling) monovalent
ions originating from a flat uniformly charged and impenetrable electrode.
These profiles are only a function of the distance from the electrode. You
may know of the Poisson-Boltzmann equation, which is the more simple
formalism of this problem that assumes point-sized ions in a
mean-electrostatic field. We can move beyond this simpler picture by adding
an excess component to the chemical potential (described by the
Carnahan-Starling equation of state) to the ideal chemical potential assume=
d
by the Poisson-Boltzmann, which I detail in equation form directly below an=
d
later on demonstrate in my Mathematica code.
Non-dimensionalized equations:
Poisson's Equation
voltage''[x] = -alpha^2 (up[x] =96 um[x])
Note: I've introduced up[x] and um[x], which are defined as:
up[x] = rho_cations[x]/rho_bulk =96 1
um[x] = rho_anions[x]/rho_bulk - 1
And up[0]=um[0]=0.
This is subject to boundary conditions:
voltage[-0.5] = surface_voltage
voltage[0] = bulk_voltage = 0
For ease of viewing, let me introduce the equation for the
non-dimensionalized excess chemical potential now:
mu_ex[Phi] = Phi(8-9*Phi + 3*Phi^2) / (1 =96 Phi)^3
Where Phi is the volume fraction of equi-sized ions and the bulk volume
fraction, Phi_bulk, is a known value:
Phi[um,up]=Phi_bulk* (2 + um[x] + up[x])
Because we are solving for an equilibrate system, we know that the
derivative of the chemical potentials for both anions and cations is equal
to zero. This is represented as:
up'[x]/(1+up[x]) + voltage'[x] +Phi_bulk *(up'[x]+um'[x])*mu_ex'[Phi] =0
um'[x]/(1+um[x]) - voltage'[x] +Phi_bulk *(up'[x]+um'[x])*mu_ex'[Phi] =0
Subject to boundary conditions:
um[0]=up[0]=0
I'm using NDSolve for this problem with the following code to solve for
voltage[x] as well as um[x] and up[x] from which I can easily obtain
concentration profiles (Mathematica code also attached):
Clear[alpha, PhiBulk, SurfaceVoltage, BulkVoltage, CPExPrime, voltage, up,
um];
alpha = 25.0; (*This can vary from 24 to 25.5*)
PhiBulk = 0.22; (*This can vary from 0.00041 to 0.22*)
SurfaceVoltage = -2.27; (*This can vary from -6 to +6*)
BulkVoltage = 0.; (*This is true for all cases*)
CPExPrime[Phi_] := (2 (4 - Phi))/(-1 + Phi)^4; (*This is the analytical
derivative of the excess chemical potential*)
Timing[
sol = NDSolve[{
voltage''[x] == -(alpha^2)*(up[x] - um[x]),
up'[x]/(1 + up[x]) + voltage'[x] + CPExPrime[PhiBulk*(2 + up[x] +
um[x])]*PhiBulk*(up'[x] + um'[x]) == 0,
um'[x]/(1 + um[x]) - voltage'[x] + CPExPrime[PhiBulk*(2 + up[x] +
um[x])]*PhiBulk*(up'[x] + um'[x]) == 0,
voltage[-1] == SurfaceVoltage,
voltage[0] == BulkVoltage,
up[0] == um[0] == 0
}, {up, um, voltage}, {x, -1, 0}][[1]]
]
My problem:
These ODEs are too stiff for Mathematica for the range of parameters I am
interested in and/or how I have formulated the problem. As I am brand new t=
o
Mathematica, I am looking for any advice on how to solve this problem. Thank you.