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.