1-D PDE with nonlinear ODE as boundary condition

Question

Recently, I am trying to solve a 1-D PDE with a nonlinear boundary condition using the function NDSolveValue. However, it seems that MMA (12) cannot solve it directly with some computational issues.

The governing equation along with associated initial and boundary conditions are where C, D, and E are constant and set as 10, 1, and 1, respectively. Note that both conditions Eqs. (3) and (4) for inner boundary conditions are required. Accordingly, my code is like

c = 10; d = 1; e = 1; sys = {(1/r)*D[r*Derivative[1, 0][f][r, t], r] == Derivative[0, 1][f][r, t] + 
          NeumannValue[c*D[g[t], t], r == 1], DirichletCondition[
         f[r, t] + (d + e*Derivative[1, 0][f][r, t])*Derivative[1, 0][f][r, t] == g[t], 
         r == 1], g[0] == 1, f[r, 0] == 0, f[5, t] == 0}; 
      {fa, ga} = NDSolveValue[sys, {f, g}, {t, 0, 1000}, MaxStepSize -> 0.00001]

However, MMA said "There are more dependent variables". Is it possible to address this issue or MMA cannot deal with nonlinear PDE?

Following the suggestion of xzczd, I rearranged the code by combining (3) and (4) and the code becomes

c = 10; d = 1; e = 1; sys = {(1/r)*D[r*Derivative[1, 0][f][r, t], r] == 
     Derivative[0, 1][f][r, t], DirichletCondition[
     f[r, t] == g[t] - (d + e*c*D[g[t], t])*c*D[g[t], t], r == 1], g[0] == 1, 
    f[r, 0] == 0, f[5, t] == 0}; 
{fa, ga} = NDSolveValue[sys, {f, g}, {t, 0, 1000}, MaxStepSize -> 0.00001]

Note that the inner boundary condition is Dirichlet type only herein. But the code cannot be calculated with errors warned by MMA.

Answer 1

NDSolve currently can't handle coupled PDE and ODE, so let's discretize the system all by ourselves, but before that, I'd like to point out the system actually has two solutions. By eliminating $\frac{∂f}{∂r}$ from $(3)$ and $(4)$ we obtain:

c = 10; d = 1; e = 1;
With[{f = f[r, t], g = g[t]}, 
 bc = {D[f, r] == c D[g, t], f - (d - e D[f, r]) D[f, r] == g} /. r -> rL;
 bcR = f == 0 /. r -> rR;
 ic = {f == 0, g == 1} /. t -> 0;
 eq = D[f, t] == 1/r D[r D[f, r], r];
 newbceq = Equal @@@ Flatten@Solve[Eliminate[bc, D[f, r] /. r -> 1], D[g, t]]]
(*
{Derivative[1][g][t] == (1/20)*(1 - Sqrt[1 - 4*f[1, t] + 4*g[t]]), 
 Derivative[1][g][t] == (1/20)*(1 + Sqrt[1 - 4*f[1, t] + 4*g[t]])}
*)

It's clear each element of newbceq leads to a solution. We choose the first one to continue (pdetoode has been used for discretization):

points = 25; domain = {rL, rR} = {1, 5};
grid = Array[# &, points, domain];
difforder = 2;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[f[r, t], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = ptoofunc[eq] // del;
odebcL = ptoofunc@{bc[[1]], newbceq[[1]]}

odebcR =(*With[{sf=1},D[#,t]+sf #&/@*)ptoofunc@bcR(*]*)
odeic = {ptoofunc[ic[[1]]] // del, ic[[2]]};
{sollst, solg} = 
 NDSolveValue[{ode, odebcL, odebcR, odeic}, {f /@ grid, g}, {t, 0, 10}(*,
  SolveDelayed->True*)]

solg // ListLinePlot

solf = rebuild[sollst, grid, 2]

Plot3D[solf[r, t], {r, rL, rR}, {t, 0, 10}, PlotRange -> All, PlotPoints -> 50]

Using newbceq[[2]] instead of newbceq[[1]] in code above we can obtain another solution: