3
$\begingroup$

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

I know its solution is

enter image description here

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
6
$\begingroup$ 
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"}]

fig1

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks Alex, How I can solve this differential equation for an Analytical solution like the one I have mentioned using DSolve or something else. $\endgroup$
    – ubaid ullah
    Commented Jul 22, 2018 at 3:35
2
$\begingroup$

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]

enter image description here

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.

Add separation of variables solution.

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

enter image description here

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

enter image description here

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]

enter image description here

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – ubaid ullah
    Commented Jul 23, 2018 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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