0
$\begingroup$

I would like to solve the following 2D+1 PDE

$\partial_t \mathcal{P}(t,\theta,\omega)=-\partial_{\theta}( \omega\mathcal{P})+\partial_{\omega}\left[(\tau_i n\cos(n\theta)-\tau_m(\theta))\mathcal{P}\right]+\partial_{\omega}\left[\left(\gamma(\theta)\omega+\gamma(\theta)\partial_{\omega}\right)\mathcal{P}\right]$

The main problem is $\tau_m(\theta)$ and $\gamma(\theta)$ are complicated functions. Is there any way to make these functions stored in advance just like $\cos(\theta),\sin(\theta)$ to make the calculation faster?

The below are $\tau_m(\theta), \gamma(\theta)$, which are just simple periodic functions.

enter image description here

Here is my original code

n = 3; 
φ = π/2;
vg0 = 5;
(*Tp1=15;Tp2=30;T=Tp2+30;*)
vg[t_] := vg[t] = vg0;
(*vg[t_]:=vg[t]=2vg0*1/(E^(k(t-Tp1))+1)-2vg0*1/(\
E^(k(t-Tp2))+1)+vg0;*)

τi = 3; Γ = 10; k = 2;

Vb = 0; µ[α_] := (-1)^(α + 1) Vb/2  ;


XTicks1 = Table[2 π*j, {j, -10, 10}];
XTicks2 = Table[π/6*j, {j, -10, 10}];
YTicks = Table[2 π*j, {j, -10, 10}];



(*//////////////////////////////////////         Karrasch poles and \
coefficients          \
/////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

Np = 25; M = 2 Np;
B = Normal[
   SparseArray[{Band[{2, 1}] -> 
      Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}], 
     Band[{1, 2}] -> 
      Table[N[1/(2 Sqrt[(2 n - 1) (2 n + 1)])], {n, 1, M - 1}]}, M]];
{bvals, bvecs} = Eigensystem[B];
Zp = Table[Abs[N[1/bvals[[2 p]]]], {p, 1, Np}];
Rp = Table[
   N[(Normalize[bvecs[[2 p]]][[1]]/(2 bvals[[2 p]]))^2], {p, 1, Np}];

(*/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////\
*)

σ0[θ_, V_, τ0_, Γ_, 
   Vg_, φ_] := σ0[θ, V, τ0, Γ, 
    Vg, φ] = 
   1/2 - I/(
     4 π) (PolyGamma[
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] - 
       PolyGamma[
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)] - 
       PolyGamma[
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)]);
τ[θ_, V_, τ0_, Γ_, 
   Vg_, φ_] := τ[θ, V, τ0, Γ, 
    Vg, φ] = -τ0 n (Cos[n θ + φ] - 
      Cos[n  θ]) σ0[θ, 
     V, τ0, Γ, Vg, φ];
σ1[θ_, V_, τ0_, Γ_, 
   Vg_, φ_] := σ1[θ, V, τ0, Γ, 
    Vg, φ] = 
   1/(8 π^2 Γ)*(PolyGamma[1, 
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) - V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 - I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)] + 
       PolyGamma[1, 
        1/2 + I/(
          2 π) (τ0 (Sin[n θ + φ] - 
               Sin[n θ]) + V/2 + Vg) + Γ/(
         4 π)]) - \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(p = 1\), \(Np\)]\(Rp[\([p]\)] \((
\*FractionBox[\(3 
\*SuperscriptBox[\((τ0 \((Sin[n\ θ + φ] - \ 
                 Sin[n\ θ])\) - V/2 + 
              Vg)\), \(2\)] \((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\) - 
\*SuperscriptBox[\((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\), \(3\)]\), 
SuperscriptBox[\((
\*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + φ] - \ 
                  Sin[n\ θ])\) - V/2 + Vg)\), \(2\)] + 
\*SuperscriptBox[\((Γ/2 + 
               Zp[\([\)\(p\)\(]\)])\), \(2\)])\), \(3\)]] + 
\*FractionBox[\(3 
\*SuperscriptBox[\((τ0 \((Sin[n\ θ + φ] - \ 
                 Sin[n\ θ])\) + V/2 + 
              Vg)\), \(2\)] \((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\) - 
\*SuperscriptBox[\((Γ/2 + 
             Zp[\([\)\(p\)\(]\)])\), \(3\)]\), 
SuperscriptBox[\((
\*SuperscriptBox[\((\ τ0 \((Sin[n\ θ + φ] - \ 
                  Sin[n\ θ])\) + V/2 + Vg)\), \(2\)] + 
\*SuperscriptBox[\((Γ/2 + 
               Zp[\([\)\(p\)\(]\)])\), \(2\)])\), \(3\)]])\)\)\);
Uprime[θ_, τ0_, φ_] := τ0 n Cos[n θ];

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

a = 1;
T = 10;
ωb = -5;
ωt = 5;
θ0 = -π/4;
τm[θ_] := τm[θ] = 
   N[τ[θ, Vb, τi, Γ, vg0, φ]];
γ[θ_] := γ[θ] = 
   N[(τi n (Cos[n θ + φ] - 
         Cos[n  θ]))^2 σ1[θ, 
      Vb, τi, Γ, vg0, φ]];
Plot[{γ[θ], τm[θ]}, {θ, -π, \
π}]
With[{u = u[t, θ, ω]}, 
  eq = D[u, 
     t] == -ω D[ 
       u, θ] + (τi n Cos[
          n θ] - τm[θ]) D[
       u, ω] + γ[θ] D[ ω u, ω] + \
γ[θ] D[ u, {ω, 2}];
  ic = u == 
     EllipticTheta[3, (θ - θ0)/2, E^(-a^2/2)]*
      E^(-ω^2/(2 a^2))/(2  π)^(3/2) a /. t -> 0];

ufun = NDSolveValue[{eq, ic, 
     u[t, -π, ω] == u[t, π, ω], 
     u[t, θ, ωb] == 0, u[t, θ, ωt] == 0}, 
    u, {t, 0, 
     T}, {θ, -π, π}, {ω, ωb, ωt}, 
    Method -> mol[35], MaxSteps -> Infinity]; // AbsoluteTiming
plots = Table[
    Plot3D[Abs[
      ufun[t, θ, ω]], {θ, -π, π}, {\
ω, ωb, ωt}, PlotRange -> All, 
     AxesLabel -> Automatic, PlotPoints -> 50, 
     BoxRatios -> {Pi, ωb, 1}, 
     ColorFunction -> "LakeColors"], {t, 0, T, 1}]; // AbsoluteTiming
ListAnimate[plots] // AbsoluteTiming
$\endgroup$
8
  • 1
    $\begingroup$ Memoization for functions with real arguments is usually rather pointless as the chances that exactly the same pops up somewhere else are rather low. Apart from that: Do you really expect us to go through this entire convoluted code and to find out which parts are particularly slow? $\endgroup$ Aug 21, 2018 at 13:07
  • $\begingroup$ My point is the function itself is very simple once we evaluate it. I am wondering if one can access the function by interpolation instead of function definition. $\endgroup$
    – Bob Lin
    Aug 21, 2018 at 13:41
  • $\begingroup$ That could be an option. We don't you just try it? $\endgroup$ Aug 21, 2018 at 13:44
  • $\begingroup$ That's my question. Is there any built-in function dealing with this? $\endgroup$
    – Bob Lin
    Aug 21, 2018 at 13:46
  • 2
    $\begingroup$ The functions $\sigma 0, \sigma 1, \tau$ are computed fairly quickly throughout the solution domain. They do not slow down the solution of the problem, but the code itself is written horribly, it needs to be rewritten from the beginning to the end. $\endgroup$ Aug 21, 2018 at 14:47

0

Your Answer

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

Browse other questions tagged or ask your own question.