Re: Stiff ODE: Modified Poisson-Boltzmann
- To: mathgroup at smc.vnet.net
- Subject: [mg118653] Re: Stiff ODE: Modified Poisson-Boltzmann
- From: Brian Giera <bgiera at gmail.com>
- Date: Fri, 6 May 2011 07:22:12 -0400 (EDT)
- References: <BANLkTimbzDTaQVfb5z7CZeEpvA6dE6Nw=w@mail.gmail.com>
Apologies, I need to update my question with the following code. Having
double checked my non-dimensionalization, I realize I am off by a factor of
1/2, which changes the range of values for alpha and one boundary condition.
With this change, I am still suffering from a stiff ODE.
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]]
]
On 05/05/11 19:56, Brian Giera wrote:
Hello,
I have included my question below. I am new to Mathematica as well as
the newsgroup, but I assume any replies will come back to me through
this email address (as opposed to me needing to check on the website).
Please let me know if I need to fix anything with my post.
Thanks,
Brian
----
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 assumed by the Poisson-Boltzmann,
which I detail in equation form directly below and later on demonstrate
in my Mathematica code.
Non-dimensionalized equations:
Poisson's Equation
oltage''[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 = 50.0; (*This can vary from 48 to 51*)
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*)
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[-0.5] == SurfaceVoltage,
voltage[0] == BulkVoltage,
up[0] == um[0] == 0
}, {up, um, voltage}, {x, -0.5, 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 to Mathematica, I am looking for any advice on how to solve this problem.
Thank you.