# Spaced time mixed partial differential equation

I am trying to reproduce solution of the following differential equation.

I know its solution is

I tried the following code.

equ = NDSolve[{D[u[\[Tau], r], {\[Tau], 1}] -
D[u[\[Tau], r], {r, 2}] - ((1/r)*D[u[\[Tau], r], {r, 1}]) == 0},
Derivative[0, 1][u][\[Tau], 0] == 0,
Derivative[0, 1][u][\[Tau], RC] == -h*u, u[0, r1] = u0,
u[0, r2] = 0, {0 <= r1 >= R}, {R <= r1 >= RC}, {r, 0, RC}, {\[Tau],
0, Infinity}]


but ended with the error "NDSolve::dsvar:$u^{0,1}[\tau,RC]=-hu$ cannot be used as a variable. I want to solve this equation keeping both $h$ and $u$ as variables and reproduced the given solution. Currently, I am out of idea how to deal with it.

Waiting for some valuable directions and suggestions.

Note: In my code a=R and b=RC.

x0 = 10^-10; xN = 2; a = 1; t0 = 1; h = .5; u0 = 1;
eq = D[u[t, x], t] - D[u[t, x], x, x] - D[u[t, x], x]/x == 0;
ic = u[0, x] == Piecewise[{{u0, x < a}}, 0];
bc = {Derivative[0, 1][u][t, x0] == 0,
Derivative[0, 1][u][t, xN] == -h*u[t, xN]};
sol = NDSolveValue[{eq, ic, bc}, u, {t, 0, t0}, {x, x0, xN}]; // Quiet
Plot3D[sol[t, x], {t, 0, t0}, {x, x0, xN}, Mesh -> None,
ColorFunction -> Hue, AxesLabel -> {"t", "x"}]


• Thanks Alex, How I can solve this differential equation for an Analytical solution like the one I have mentioned using DSolve or something else. Commented Jul 22, 2018 at 3:35

Your bc's set up nicely, so we can solve this pde as a function of u[t,r]. The problem is the singular term u'[r]/r at r = 0, but since you have specified that u'[0] == 0 we have 0/0 at r = 0 so we can use L'Hospital's rule to get the limit there. Taking derivatives of the numerator and denominator we see that u'[r]/r -> u''[r]/1 for r = 0. Define the pde in Piecewise form as:

pde = D[u[t, r], {r, 2}] + Piecewise[{{D[u[t, r], {r, 2}], r == 0}, {D[u[t, r], r]/r, r > 0}}] == D[u[t, r], t]


Now NDSolve has no problem at r = 0.

the bc's are:

bc1 = Derivative[0, 1][u][t, 0] == 0
bc2 = Derivative[0, 1][u][t, RC] == -h*u[t, RC]


the initial condition

ic = u[0, r] == Piecewise[{{u0, r < R}}, 0]


Plug in some values so NDSolve can actually do something.

R = 1;
RC = 3 R;
h = 1;
u0 = 1;

sol = NDSolve[{pde, bc1, bc2, ic}, u[t, r], {t, 0, 10}, {r, 0, RC}];

u[t_, r_] = u[t, r] /. %[[1]];


Plot the solution.

tb = Table[Plot[Evaluate[u[t, r]], {r, 0, RC}, PlotRange -> {0, 1}], {t, 0, 5, .1}];

ListAnimate[tb]


I don't know if there is a way that DSolve can solve the problem analytically directly, but it should be relatively easy through separation of variables.

At first look for separable solutions of the type u = R[r]*T[t]

Clear["Global*"]

pde = D[u[r, t], {r, 2}] + (1/r)*D[u[r, t], r] == D[u[r, t], t]

u[r_,t_]=R[r]T[t]

pde = Expand[pde[[1]]/u[r, t] == pde[[2]]/u[r, t]]

R''[r]/R[r]+R''[r]/(r R[r])==T''[t]/T[t]


Left side is dependent only or r and right side is dependent only on t. So each side must be equal to the same constant. That constant can be 0, >0, or <0, each case giving a different set of solutions. In this case a negative constant is the only one we need.

$Assumptions= c>0 teq=pde[[2]]==-c^2 DSolve[teq,T[t],t]//Flatten; T1=(T[t]/.%/.C[1]->1); req= pde[[1]]==-c^2; DSolve[req,R[r],r]//Flatten; R1=(R[r]/.%/.{C[1]->c1,C[2]->c2});  The solution then becomes: un[r_,t_]=R1 T1 (*E^(-c^2 t) (c1 BesselJ[0,c r]+c2 BesselY[0,c r])*)  Recognizing that we will have an infinite series of un's. The solution must be finite at origin, so c2=0  Apply the boundary condition bc=(D[un[r, t], r] /. r -> b) + h*un[b, t] == 0//Collect[#, {c1, E^(-t c^2)}]& (*c1 E^(-c^2 t) (h BesselJ[0,b c]-c BesselJ[1,b c])==0*)  This equation sets the values of c. There are an infinite number of them and they are the roots of: ceq=(h BesselJ[0,b c]-c BesselJ[1,b c])==0  Solve the above for BesselJ[1,b c] so we can substitute later. csol=Solve[ceq,BesselJ[1,b c]]//Flatten;  Use the ic for u[r,0] to get c1. We will use the orthogonality of BesselJ0. The ic is easiest( for me) to define in terms of UnitStep functions. ic=un[r,0]==u0(UnitStep[r]-UnitStep[r-a]) (*c1 BesselJ[0,c r]==u0 (UnitStep[r]-UnitStep[r-a])*)  Multiply both sides of the above by r BesselJ0[r c] and Integrate over r from a to b $Assumptions=\$Assumptions&&a>0&&b>a

Integrate[ic[[1]]*r*BesselJ[0, c*r], {r, 0, b}] == Integrate[ic[[2]]*r*BesselJ[0, c*r], {r, 0, b}];


Apply the condition for c

%/.csol//Simplify

Solve[%,c1]//Flatten;

c1=c1/.%


We end up with

un[r,t]/.c->c[[n]]

(2 a u0 c[[n]] E^(-t c[[n]]^2) BesselJ[1,a c[[n]]] BesselJ[0,r c[[n]]])/(b^2 (c[[n]]^2+h^2) BesselJ[0,b c[[n]]]^2)


The actual answer as an infinite sum over the c's.

u[r_, t_] := ((2*a*u0)/b^2)*Sum[(c[[n]]*BesselJ[1, a*c[[n]]]*BesselJ[0, r*c[[n]]])/
(E^(t*c[[n]]^2)*((c[[n]]^2 + h^2)*BesselJ[0, b*c[[n]]]^2)), {n, 1, Infinity}]


Plug in some numbers so we can plot.

a=1;
b=3 a;
h=1;
u0=1;


Solve for the c's from

F1[c_]=ceq[[1]]
(*BesselJ[0,3 c]-c BesselJ[1,3 c]*)


I really tried to use NSolve in the following form to solve for the c's, because it has worked before on other problems:

c=d/.NSolve[{F1[d],d>=0&&d<=100},d];


I got an array of most of the values, but some were missing. Increasing the WorkingPrecision up to 200 added some but still did not get all of them. So using FindRoot manually, first plot the first two values to get a starting point.

Clear[c];
Plot[F1[c],{c,0,2}]


Find 200 values

n = 200;
c = Table[0, {i, n}];
aa = FindRoot[F1[\[Beta]] == 0, {\[Beta], 0.6}, WorkingPrecision -> 20, AccuracyGoal -> 15];
c[[1]] = \[Beta] /. aa;
aa = FindRoot[F1[\[Beta]] == 0, {\[Beta], 1.5}, WorkingPrecision -> 20, AccuracyGoal -> 15];
c[[2]] = \[Beta] /. aa;
For[i = 3, i <= n, aa = FindRoot[F1[\[Beta]] == 0, {\[Beta], 2*c[[i - 1]] - c[[i - 2]]}, WorkingPrecision -> 20, AccuracyGoal -> 15];
c[[i]] = \[Beta] /. aa; i++]
Clear[n]

u[r_, t_, m_] := ((2*a*u0)/b^2)*Sum[(c[[n]]*BesselJ[1, a*c[[n]]]*BesselJ[0, r*c[[n]]])/
(E^(t*c[[n]]^2)*((c[[n]]^2 + h^2)*BesselJ[0, b*c[[n]]]^2)), {n, 1, m}]


We have slow convergence for u at t = 0 since Exp[-c^2 t] does not contribute. Check that we satisfy the initial condition with the 200 terms.

Plot[Evaluate[u[r,0,200]],{r,0,b}]


Matches the ic within numerical tolerance. For t > 0 the series converges with many fewer terms.

tb1=Table[Plot[Evaluate[u[r,t,20]],{r,0,b},PlotRange->{0,1}],{t,0,5,.1}];
ListAnimate[tb1]


It's pretty much the same as the NDSolve` solution, but the analytical solution was certainly not as quick.

• Thanks Bil Watts, for this nice valuable comment. I am trying to solve this on paper pen if failed then I will rise another question as u suggested. Commented Jul 23, 2018 at 11:23