DSolve[] is unable to compute the symbolic solution of a non-linear BVP with double Neumann BCs

Question

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]

Answer 1

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.

Addendum

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]"}]