I am looking for a symbolic closed-form solution to the following BVP.

$$\frac{\partial ^2}{\partial x}\phi(x) = k j(x)$$

where $k$ is a real-valued scalar constant, and $j(x)$ is a strongly non-linear function given by $$j(x) = j_0 (e^{\frac{(1-\alpha)F \eta(x)}{RT}} - e^{\frac{-\alpha F \eta(x)}{RT}})$$

wherein, $j_0$, $\alpha$, $R$ & $T$ are scalar constants, and

$$\eta(x) = \phi(x) - a - b$$

where $a$ and $b$ are two other constants.

The boundary conditions are $\phi'(0) = 0$ & $\phi'(1) = 1$ in the unit domain $(0,1)$

I tried solving this using DSolve, but was unsuccessful. Then, I tried solving for a linearised version of $j(x)$ using a Taylor expansion about an operating point $\phi^*$. This code works and produces a symbolic answer. I do realise that non-linear problems are hard, and do not have a general solution. But in this case, I do know that we have a solution (from theory, and also was able to compute it numerically with a relative error/tolerance of $10^{-7}$).

Furthermore, even though we have a double Neumann BC for this elliptic problem, we have a unique solution due to the fact that the source term is implicit, i.e. it contains the field variable $\phi(x)$ being solved for.

Any help regarding solving this problem in Mathematica is highly appreciated. The code I tried for the linearised version is given below

Eta =  Phi[x] - a - b
linearisedBV = Collect[Normal[Series[(Exp[((1 - alpha) F Eta)/(R T)] -  Exp[ (-alpha F Eta)/(R T)]), {Phi[x], phistar, 1}]], Phi[x]]
eqn = Phi''[x] == k*j0*linearisedBV
soln = DSolve[{eqn, Phi'[0] == 0, Phi'[1] == 1}, Phi[x], x]
  • $\begingroup$ soln = DSolve[{eqn, Phi[0] == 0, Phi'[1] == 1}, Phi[x], x] also, you're missing a curly bracket before Phi[x] in linearisedBV $\endgroup$
    – Feyre
    Aug 22 '16 at 12:28
  • $\begingroup$ Yup.Sorry, I have the correct syntax in my mathematica file and it runs fine. They are just typos that crept in when I copied and pasted onto this forum. The typos in the code are now fixed. $\endgroup$ Aug 22 '16 at 13:01
  • $\begingroup$ The last line of the linearized calculation should contain Phi'[0] == 0 instead of Phi[0] == 0. $\endgroup$
    – bbgodfrey
    Aug 23 '16 at 5:15

The nonlinear problem in the question,

Eta =  Phi[x] - a - b;
eqn = Phi''[x] == k*j0*(Exp[((1 - alpha) F Eta)/(R T)] -  Exp[ (-alpha F Eta)/(R T)]);
soln = DSolve[{eqn, Phi'[0] == 0, Phi'[1] == 1}, Phi[x], x]

yields an empty solution set with the explanatory warning

DSolve::bvimp: General solution contains implicit solutions. In the boundary value problem, these solutions will be ignored, so some of the solutions will be lost.

DSolve, or more precisely Solve, which it calls, cannot invert the boundary conditions. Sometimes,

SetOptions[Solve, Method -> Reduce]

overcomes this difficulty, but not here. Instead, apply DSolve without boundary conditions, which yields,

s = soln = DSolve[eqn, Phi[x], x]

(* Solve[
     Integrate[1/Sqrt[(2*E^((F*(a*(-1 + alpha) + (-1 + alpha)*b - alpha*K[1]))/(R*T))*
     ((-1 + alpha)*E^(((a + b)*F)/(R*T)) - alpha*E^((F*K[1])/(R*T)))*j0*k*R*T)/
     ((-1 + alpha)*alpha*F) + C[1]], {K[1], 1, Phi[x]}]^2 == (x + C[2])^2, Phi[x]] *)

That neither Solve nor Reduce can make headway with this equation is, perhaps, not surprising. Nonetheless, some progress can be made. First, extract the equation from within s and take the square root of both sides.

ss = s[[1]] /. Power[z_, 2] -> z

(* Integrate[1/Sqrt[(2*E^((F*(a*(-1 + alpha) + (-1 + alpha)*b - alpha*K[1]))/(R*T))*
   ((-1 + alpha)*E^(((a + b)*F)/(R*T)) - alpha*E^((F*K[1])/(R*T)))*j0*k*R*T)/
   ((-1 + alpha)*alpha*F) + C[1]], {K[1], 1, Phi[x]}] == x + C[2] *)

(Edit: The following was rewritten to properly address boundary conditions.) The boundary conditions now can be applied to determine the constant C[1] as follows. To obtain Phi'[x], differentiate ss, multiply the resulting equation by the denominator of the left side of the equation, and square both sides.

sss = (Numerator[D[ss, x][[1]]])^2 == (Denominator[D[ss, x][[1]]])^2

(* Derivative[1][Phi][x]^2 == 
   (2*E^((F*(a*(-1 + alpha) + (-1 + alpha)*b - alpha*Phi[x]))/(R*T))*
   ((-1 + alpha)*E^(((a + b)*F)/(R*T)) - alpha*E^((F*Phi[x])/(R*T)))*j0*k*R*T)/
   ((-1 + alpha)*alpha*F) + C[1 *)

Note that this result also can be obtained by multiplying eqn by Phi'[x] and integrating. Now C[1] can be computed trivially by

First@Solve[sss /. x -> 0 /. Phi'[0] -> 0, C[1]]

(* {C[1] -> (-2*E^((F*(a*(-1 + alpha) + (-1 + alpha)*b - alpha*Phi[0]))/(R*T))*
    (-E^(((a + b)*F)/(R*T)) + alpha*E^(((a + b)*F)/(R*T)) - alpha*E^((F*Phi[0])/(R*T)))*j0*
    k*R*T)/((-1 + alpha)*alpha*F)} *)

as well as from First@Solve[sss /. x -> 1 /. Phi'[1] -> 1, C[1]]. Consistency of these two results requires

1 == (sss[[2]] /. x -> 1) - (sss[[2]] /. x -> 0)

(* 1 == (-2*E^((F*(a*(-1 + alpha) + (-1 + alpha)*b - alpha*Phi[0]))/(R*T))*
   ((-1 + alpha)*E^(((a + b)*F)/(R*T)) - alpha*E^((F*Phi[0])/(R*T)))*j0*k*R*T)/
   ((-1 + alpha)*alpha*F) + 
   (2*E^((F*(a*(-1 + alpha) + (-1 + alpha)*b - alpha*Phi[1]))/(R*T))*
   ((-1 + alpha)*E^(((a + b)*F)/(R*T)) - alpha*E^((F*Phi[1])/(R*T)))*j0*k*R*T)/
   ((-1 + alpha)*alpha*F) *)

A solution to the nonlinear ODE eqn therefore exists only if this consistency condition can be satisfied. Adjusting C[2] may well be sufficient to satisfy the consistency condition, but determining the necessary value for C[2] appears to be a numerical calculation.


A little more progress can be made for the special case, alpha == 1/2.

eqn = Phi''[x] == 2 k*j0*Sinh[Eta F/(2 R T)];
soln = Phi[x] /. DSolve[{eqn}, Phi[x], x] /. 
    JacobiAmplitude[z1_, z2_] :> JacobiAmplitude[Simplify[z1], z2]

(* {(a F + b F + 4 I R T JacobiAmplitude[(I Sqrt[F] Sqrt[8 j0 k R T + F C[1]] (x + C[2]))/
       (4 R T), (16 j0 k R T)/(8 j0 k R T + F C[1])])/F, 
    {(a F + b F - 4 I R T JacobiAmplitude[(I Sqrt[F] Sqrt[8 j0 k R T + F C[1]] (x + C[2]))/
       (4 R T), (16 j0 k R T)/(8 j0 k R T + F C[1])])/F} *)

Unfortunately, the resulting boundary conditions cannot be solved symbolically by Mathematica for the two constants of integration.

Solve[{0 == D[soln[[1]], x] /. x -> 0, 1 == D[soln[[1]], x] /. x -> 1}, {C[1], C[2]}]

returns unevaluated with the message

Solve::nsmet: This system cannot be solved with the methods available to Solve.

Reduce does no better.

Second Addendum

The OP indicates in a comment below that only Phi[0] is needed. I know of no way to obtain even that symbolically. However, it can be obtained numerically without difficulty. For instance, with a + b == 0, k j0 == 1, and F/(R T) == 1,

eqn = Phi''[x] == Exp[(1 - alpha) Phi[x]] - Exp[ -alpha Phi[x]];
ParametricNDSolveValue[{eqn, Phi'[0] == 0, Phi'[1] == 1}, Phi[0], x, {alpha}];
Plot[%[alpha], {alpha, 0, 1}, AxesLabel -> {alpha, "Phi[0]"}]

enter image description here

  • $\begingroup$ I have understood your derivation and code thus far, including the calculation of C[1] and the corresponding consistency condition. However, I am still unable to understand how one might proceed further with this calculation. How do we link this back to C[2] (as per the last sentence of your post) ? If it will help any further, I am NOT interested in solution over the entire domain. All that I need is the value of Phi[0]. What are the potential next steps you'd recommend to solve this problem ? $\endgroup$ Aug 25 '16 at 12:05
  • $\begingroup$ @Krishna I believe that progress can be made only if numerical values are assigned to the three constants alpha, a + b, and F/(R T). $\endgroup$
    – bbgodfrey
    Aug 25 '16 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.