# How to accelerate NDSolve for PDE with user defined functions?

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.

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

• 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? Aug 21, 2018 at 13:07
• 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. Aug 21, 2018 at 13:41
• That could be an option. We don't you just try it? Aug 21, 2018 at 13:44
• That's my question. Is there any built-in function dealing with this? Aug 21, 2018 at 13:46
• 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. Aug 21, 2018 at 14:47