# 1-D PDE with nonlinear ODE as boundary condition

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.

• what is the equation for g? If there is none, then g can not be a dependent variable. Also, note that you are using the finite element method and not a finite difference method. See here and the paragraph above it. I am just mentioning this because you added the FDM tag. Jun 4 '20 at 7:33
• your code does not run. gives syntax error. But do you really need all this other stuff , like timing and monitor and all this in order to show what you think is an issue with Boundary conditions? Why not make a MWE that just shows the essential part of the problem. Anything else you can remove as not needed. Jun 4 '20 at 7:36
• Seems that you've copied the code in wrong way. (You should not select "Copy As -> Plain Text". Things like Derivative[……] won't be copied correctly in this way. ) Please Ctrl+Shift+I to convert the code to input form, and Ctrl+C to copy it, for more info check this post: mathematica.meta.stackexchange.com/a/1585/1871 @user21 I think we can eliminate $g$ using $(3)$ and $(4)$. Jun 4 '20 at 7:45
• Thank you. I changed the format of the code. Jun 4 '20 at 7:56
• @xzczd if eliminate g using both (3) and (4), how to preserve the initial condition in Eq. (6)? Jun 4 '20 at 7:56

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,
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:

• What's an amazing work! I'll try this. Jun 4 '20 at 9:57