1
$\begingroup$

I have the following equation $$\frac{\partial\,P(\theta,t)}{\partial t} = \alpha\cos\theta\,P(\theta,t) + \beta \frac{\partial^2\,P(\theta,t)}{\partial \theta^2} $$ subject to following conditions$$P(\theta, 0)=\delta(\theta-\pi/4), P(\theta,t)=P(\theta+2\pi,t)~.$$

I tried solving this numerically in Mathematica:

alpha = 0.02
beta = 0.02
fokkerPlanck = {D[p[x, t], t]==alpha*Cos[x]*p[x,t] + beta*D[p[x,t],{x,2}], p[x,0] == DiracDelta[Pi/4], p[2Pi, t] == p[0, t]};
sol = Flatten@NDSolve[fokkerPlanck, p, {x, 0, 2Pi},{t, 0, 100}]

If I run this, I get null solution, so something is weird. Could anyone let me know what is potentially wrong here?

Improve this question
$\endgroup$
4
  • $\begingroup$ From the first glance, you need to fix also the interval in x and the corresponding boundary conditions. Am I wrong? $\endgroup$
    – Alexei Boulbitch
    Commented Jun 17, 2021 at 18:48
  • $\begingroup$ @AlexeiBoulbitch You are right that I need to add interval in x. I was doing that in my code, but while typing it here, I missed it. This is now fixed. The boundary conditions are already included in fokkerPlanck though. $\endgroup$
    – titanium
    Commented Jun 17, 2021 at 18:54
  • $\begingroup$ There is an example Fokker-Panck equation in the reference. Maybe that's useful. See also this $\endgroup$
    – user21
    Commented Jun 18, 2021 at 5:29
  • $\begingroup$ @titanium Initial condition is not periodic, while boundary condition is periodic. What kind of solution do you expect? $\endgroup$
    – Alex Trounev
    Commented Jun 18, 2021 at 11:09

2 Answers 2

Reset to default
2
$\begingroup$

It seems the problem is the the DiracDelta boundary condition.

If I replace that with a highly peaked Gaussian, p[x, 0] == Exp[-(x - \[Pi]/4)^2/(2*0.001)^2], I get a finite solution.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Even if I do so, plot of p[x, 100] looks really weird; its amplitude is way too small. $\endgroup$
    – titanium
    Commented Jun 17, 2021 at 22:18
  • $\begingroup$ Well it's a diffusion equation, the width is supposed to get larger. Also your equation does not have a $P$ on the $\alpha$ term, but your Mathematica code does? $\endgroup$
    – SuperCiocia
    Commented Jun 17, 2021 at 22:26
  • $\begingroup$ Thanks for catching that. I corrected the problem. My point was after running the code for long enough, the distribution of p[x,100] should be a Gaussian distribution which integrates to 1 in the range of 0, 2Pi. This is because p is the probability density. At the moment I cannot recover this. $\endgroup$
    – titanium
    Commented Jun 17, 2021 at 22:29
  • $\begingroup$ @titanium Gaussian distribution defined on $-\infty <x<\infty$. It is not periodic function on $0\le x\le 2 \pi$ $\endgroup$
    – Alex Trounev
    Commented Jun 18, 2021 at 11:59
1
$\begingroup$

First step is come out from singularity at $t=0$. For this we can use analytical solution in the form

sol1 = 
 DSolve[{D[p[x, t], t] == 
    1/50*Cos[x]*p[x, t] + 1/50*D[p[x, t], {x, 2}], 
   p[x, 0] == DiracDelta[(x - Pi/4)]}, p[x, t], {x, t}]

Out[]= {{p[x, t] -> 
   ConditionalExpression[(
    5 E^(-((25 (\[Pi] - 4 x)^2)/(32 t)) + 1/50 t Cos[x]))/(
    Sqrt[2 \[Pi]] Sqrt[t]), Re[t] > 0]}}

This solution is not periodic on $0\le x \le 2 \pi$, but we can make it periodic for $t>t_0$ using FEM as follows

Needs["NDSolve`FEM`"]
alpha = 0.02;
beta = 0.02; reg = Rectangle[{0, 1/2}, {2 Pi, 100}]; mesh = 
 ToElementMesh[reg, 
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    area > 0.001 (0.05 + 3 Norm[Mean[vertices]])]]
fokkerPlanck = {D[p[x, t], t] == 
    alpha*Cos[x]*p[x, t] + beta*D[p[x, t], {x, 2}], 
   DirichletCondition[
    p[x, t] == (
     5 E^(-((25 (\[Pi] - 4 x)^2)/(32 t)) + 1/50 t Cos[x]))/(
     Sqrt[2 \[Pi]] Sqrt[t]), t == 1/2 && 10^-2 <= x <= 2 Pi - 10^-2], 
   PeriodicBoundaryCondition[p[x, t], x == 2 \[Pi], 
    Function[x, x - 2 \[Pi]]]};

sol = NDSolve[fokkerPlanck, p, Element[{x, t}, mesh]]

There is a message about stability of solution

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

Visualization Figure 1

Improve this answer
$\endgroup$
4
  • $\begingroup$ If I execute the DSolve part for obtaining analytical solution, it returns me the code I entered, not the analytical result. Any idea on why that might be the case? I cleared all the cache as well, but that did not help. $\endgroup$
    – titanium
    Commented Jun 18, 2021 at 16:50
  • $\begingroup$ @titanium What version do you run? $\endgroup$
    – Alex Trounev
    Commented Jun 18, 2021 at 21:09
  • $\begingroup$ It is version 12.0.0.0. $\endgroup$
    – titanium
    Commented Jun 18, 2021 at 22:13
  • $\begingroup$ @titanium Ok! My $Version is 12.3.0 for Microsoft Windows (64-bit) (May 10, 2021). But DSolve has been updated in 12. Could you show code you execute? $\endgroup$
    – Alex Trounev
    Commented Jun 19, 2021 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.