1
$\begingroup$

\begin{align} u_t&=v*(u_{rr}+r^{-1}u_r), \quad 0\leq r <a , t>0\\ u(r,0)&=f(r)\\ u(a,t)&=0 ,\quad t \geq0 \end{align}

We know that $f(r)=J_0(x_{01}r)$, where $x_{01}$ is the first positive root of the Bessel Function and $a=v=1$

Now we solve the problem with the method separation of variables. Let $$u(r)=R(r)T(t).$$ Then, $$R(r)T'(t)=R''(r)T(t)+r^{-1} R'(r)T(t)$$ We divide both sides with $R(r)T(t)$. Thus, $$\frac{T'(t)}{T(t)}=\frac{R''(r)}{R(r)}+r^{-1}\frac{R'(r)}{R(r)}=-\lambda$$. So we have the following two equations $$r\cdot R''(r)+R'(r)+\lambda\cdot r R(r)=0, \quad 0<r<1 \tag 1$$ It is Sturm Liouville problem, from which we can conclude the weight function is $w=r$. We change variables $x=k \cdot r$ then we obtain $$x\bar{R''}+bar{R'}+x\bar{R}=0 \tag2$$. The equation $(2)$ has as solution the Bessel function $\bar{R}=J_0(x)=J_0(kr)$.

The boundary condition $u(a,t)=0, t \geq 0$ become $R(a)=0$, so $J_0(ka)=0$. That means $k \cdot a$ must be the root of the Bessel function with zero order.

$$k_m \cdot a= x_{0m}, \quad m=1,2,...$$.

Now $W'(t)+\lambda W(t)=0$ has general solution $$W(t)=c_1 e^{-\lambda \cdot t} \tag3$$ Then we have $$W(t)=c_1 e^{-k_n^2 \cdot t}$$ Now from superposition, the general solution to the problem will be $$u(r,t)=\sum_{n=1}^{\infty} d_m \cdot e^{-k_n^2 \cdot t} \cdot J_0(k_m r)$$ and the coefficient will be $$d_m=\frac{2 \int_{0}^{a} r \cdot f(r) \cdot J_0(k_m r) dr}{a^2 J_1^2(x_m )}, where \quad k_m=\frac{x_{0m}}{a}$$ Now I am trying to write the Mathematica code

a = 1;
x0[m_] := N[BesselJZero[0, m]]
(*d definition of the function f(r)*)

f[r_] := BesselJ[0, x0[1]*r]
k[m_] := x0[m]/a
(* Computation of the Fourier coeffiecients*)
d[m_] := (2 \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(a\)]\(r*
    BesselJ[0, k[m]*r]\ *f[r] \[DifferentialD]r\)\))/(
 a^2*BesselJ[1, x0[m]]^2)
(*Construction of the exact solution*)

u[r_, t_, m_] := d[m]*E^(-k[m]^2*t)*BesselJ[0, k[m]*r]
uapprox[r_, t_] = \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(n = 1\), \(10\)]\(u[r, t, n]\)\);
(*Graph of the exact solution*)
Manipulate[
 ParametricPlot3D[{r*Cos[\[Theta]], r*Sin[\[Theta]], 
   uapprox[r, t]}, {r, 0, 1}, {\[Theta], 0, 2 Pi}, 
  PlotRange -> {-1.2, 1.2}, BoundaryStyle -> Directive[Red, Thick], 
  ColorFunction -> "SolarColors", Mesh -> True], {t, 0, 10, 0.0001}]

graphs = 
 Table[ParametricPlot3D[{r*Cos[\[Theta]], r*Sin[\[Theta]], 
    uapprox[r, t]}, {r, 0, 1}, {\[Theta], 0, 2 Pi}, 
   PlotRange -> {-1.2, 1.2}, BoundaryStyle -> Directive[Red, Thick], 
   ColorFunction -> "SolarColors"], {t, 0, 2.5, 2.5/8}]

Animate[ParametricPlot3D[{r*Cos[\[Theta]], r*Sin[\[Theta]], 
   uapprox[r, t]}, {r, 0, 1}, {\[Theta], 0, 2 Pi}, 
  PlotRange -> {-1.2, 1.2}, BoundaryStyle -> Directive[Red, Thick], 
  ColorFunction -> "SolarColors", Mesh -> True], {t, 0, 10, 0.0001}]

When I try to run ParametricPlot3D gives me. Could anyone explain to me, please, why I received that answer? enter image description here

$\endgroup$
2
  • $\begingroup$ We can derive exact solution using DSolve[{D[u[t, r], t] == D[u[t, r], r, r] + D[u[t, r], r]/r, u[0, r] == BesselJ[0, BesselJZero[0, 1] r], u[t, 1] == 0}, u, {t, r}] $\endgroup$ Apr 13 at 2:58
  • $\begingroup$ @AlexTrounev and then how could I write with a proper code? Sorry but i am confused $\endgroup$ Apr 13 at 15:12

1 Answer 1

1
$\begingroup$

Due to the orthogonality of the Bessel functions, the numerator of $d_m$ is zero except for $m=1$. So $u$ is just $u(r,t)=d_1\exp\left(-x_{01}/a^2\right)J_0(x_{01}r/a)$. This is why Mathematica, having defined x0[m_] := N[BesselJZero[0, m]] numerically, gives you the precision error. If you make everything symbolic, the error disappears but now that we know simpler formula, you can just use

a=1;
b=N@BesselJZero[0,1];
u[r_,t_]:=E^(-b/a^2*t)*BesselJ[0,b*r/a];
$\endgroup$
1
  • $\begingroup$ So I have to write the following code to find the solution a = 1; b = N@BesselJZero[0, 1]; u[r_, t_] := E^(-b/a^2*t)*BesselJ[0, b*r/a]; Manipulate[ ParametricPlot3D[{r*Cos[\[Theta]], r*Sin[\[Theta]], u[r, t]}, {r, 0, 1}, {\[Theta], 0, 2 Pi}, PlotRange -> {-1.2, 1.2}, BoundaryStyle -> Directive[Red, Thick], ColorFunction -> "SolarColors", Mesh -> True], {t, 0, 10, 0.0001}] $\endgroup$ Apr 13 at 15:18

Your Answer

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

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