# How to plot a function defined in considerable time with Plot3D

I'm trying to plot a function as follows:

In[1]:=Psi1[l1_,m1_,l2_,m2_]:=Sqrt[(2*l1+1)/(4*Pi)]*LegendreP[l1,m1,Cos[theta1]]*Sqrt[(2*l2+1)/(4*Pi)]*LegendreP[l2,m2,Cos[theta2]]-Sqrt[(2*l1+1)/(4*Pi)]*LegendreP[l1,m1,Cos[theta2]]*Sqrt[(2*l2+1)/(4*Pi)]*LegendreP[l2,m2,Cos[theta1]]

In[2]:=Den1[l1_,m1_,l2_,m2_]:=1/2*Psi1[l1,m1,l2,m2]^2

In[3]:=Fun1[l1_,m1_,l2_,m2_]:=Integrate[Den1[l1,m1,l2,m2]*Sin[theta1]*Sin[theta2],{theta1,x1,Pi},{phi1,0,2*Pi},{theta2,x2,Pi},{phi2,0,2*Pi}]

In[4]:=Plot3D[Fun1[2,0,1,0],{x1,0,Pi},{x2,0,Pi},PlotRange->All]


And:

In[5]:=CumDen1[l1_,m1_,l2_,m2_]:=Plot3D[Fun1[l1,m1,l2,m2],{x1,0,Pi},{x2,0,Pi},PlotRange->All]

In[6]:=CumDen1[2,0,1,0]


But I've noticed that Mathematica takes more than 60 minutes to plot the function that if it does it in parts:

In[7]:=Fun1[2,0,1,0]

In[8]:=Plot3D[%,{x1,0,Pi},{x2,0,Pi},PlotRange->All]


With this form, Mathematica graphs the function in less than 60 seconds.

Of course, everything gets too complicated when I want to plot the function for several values of l1 and l2 as follows:

In[9]:=Do[If[l1=!=l2,Print[{l1,l2,CumDen1[l1,m1,l2,m2]}]],{l1,0,5},{l2,0,5},{m1,0,0},{m2,0,0}]


Is there any way that Mathematica can plot these functions in a considerable time?

• In this case, you might wish to consider memoizing the integration so that it does not have to be done for every point in your plot: Fun1[l1_, m1_, l2_, m2_] := Fun1[l1, m1, l2, m2] = (* stuff *) – J. M.'s ennui Aug 24 '17 at 7:22
• Simply putting an Evaluate around the function in the plot should help a lot. – mmeent Aug 24 '17 at 7:27
• In your def of Fun1 you're integrating w.r.t. phi1 and phi2. But, as far as I can tell, they don't appear in the integrand. Also, for clarity if nothing else, I would probably write Psi1 and Den1 etc. explicitly as functions of the terms you later integrate over (ie. theta1, theta2, x1, x2, etc). – aardvark2012 Aug 24 '17 at 10:46

# TL;DR;

Define

int[a_ + b_] = int[a] + int[b]
int[a_] = int[int[a, 1], 2]
int[f_ a_, i_] := f int[a, i] /; FreeQ[f, Indexed[_, i]]
int[f_, i_] := 2 Pi (Pi - Indexed[x, i]) f /; FreeQ[f, Indexed[_, i]]
int[exp[a_?(FreeQ@Indexed[_, i_]) Indexed[\[Theta], i_]], i_] =
Integrate[
Exp[a Indexed[\[Theta], i]], {Indexed[\[Theta], i], Indexed[x, i],
Pi}, {Indexed[\[Phi], i], 0, 2*Pi}]


and

Psi1[l1_, m1_, l2_, m2_] :=
Sqrt[(2*l1 + 1)/(4*Pi)]*LegendreP[l1, m1, Cos[Indexed[\[Theta], 1]]]*
Sqrt[(2*l2 + 1)/(4*Pi)]*
LegendreP[l2, m2, Cos[Indexed[\[Theta], 2]]] -
Sqrt[(2*l1 + 1)/(4*Pi)]*
LegendreP[l1, m1, Cos[Indexed[\[Theta], 2]]]*
Sqrt[(2*l2 + 1)/(4*Pi)]*
LegendreP[l2, m2, Cos[Indexed[\[Theta], 1]]]

Den1[l1_, m1_, l2_, m2_] := 1/2*Psi1[l1, m1, l2, m2]^2

Fun1[l1_, m1_, l2_, m2_] :=
Fun1[l1, m1, l2, m2] =
Simplify@ExpToTrig[
int[Expand@
TrigToExp[
Den1[l1, m2, l2, m2]*Sin[Indexed[\[Theta], 1]]*
Sin[Indexed[\[Theta], 2]] /.
Sqrt[1 - Cos[Indexed[\[Theta], i_]]^2] :>
Sin[Indexed[\[Theta], i]]] /.
E^(a_ + b_: 0) :> exp[a] exp[b]] /. exp -> Exp]

Fun1P[l1_, m1_, l2_, m2_] :=
Fun1[l1, m1, l2, m2] /. {Indexed[x, 1] -> x1, Indexed[x, 2] -> x2}


Now plot with:

Plot3D[Fun1P[2, 0, 1, 0], {x1, 0, Pi}, {x2, 0, Pi}, PlotRange -> All]


# Explanation

Apart from memoization/Evaluate, we can do a little bit more to speed up this specific function:

First, we look at some typical outputs from Den1 - you'll notice that the output only consists of products of trigonometric functions and occasionally $\sqrt{1-\cos^2\theta_i}$ (which we can simply replace with $\sin\theta_i$ since $0<\theta_i<\pi$). Now we can make our own integration function that is specialized for this type of functions:

To make things a bit easier later, we redefine Psi1 to use Indexed expressions instead of hardcoded names:

Psi1[l1_, m1_, l2_, m2_] :=
Sqrt[(2*l1 + 1)/(4*Pi)]*LegendreP[l1, m1, Cos[Indexed[\[Theta], 1]]]*
Sqrt[(2*l2 + 1)/(4*Pi)]*
LegendreP[l2, m2, Cos[Indexed[\[Theta], 2]]] -
Sqrt[(2*l1 + 1)/(4*Pi)]*
LegendreP[l1, m1, Cos[Indexed[\[Theta], 2]]]*
Sqrt[(2*l2 + 1)/(4*Pi)]*
LegendreP[l2, m2, Cos[Indexed[\[Theta], 1]]]


Next, we define our specialized integral, int. First, we make it transparent to addition:

int[a_ + b_] = int[a] + int[b]


Next, we separate integrals into ones over $\theta_1$ and $\theta_2$:

int[a_] = int[int[a, 1], 2]


where we take int[f,i] to mean "integrate f over $\theta_1,\phi_1$". Now, we need to pull out constant factors:

int[f_ a_, i_] := f int[a, i] /; FreeQ[f, Indexed[_, i]]


The Condition (/; for short) ensures that the factor we want to pull out does indeed not contain any variable indexed with i. Integrals over constants are handled by

int[f_, i_] := 2 Pi (Pi - Indexed[x, i]) f /; FreeQ[f, Indexed[_, i]]


where the Condition ensures once again that we really have a constant. To check how far we got, we can apply it to e.g. Den1[2, 1, 1, 0]*Sin[Indexed[\[Theta], 1]]*Sin[Indexed[\[Theta], 2]]:

int[Expand@
TrigToExp[
Den1[2, 1, 1, 0]*Sin[Indexed[\[Theta], 1]]*
Sin[Indexed[\[Theta], 2]]]]


We apply Expand and TrigToExp to get terms that are as simple as possible. There are still some square roots left. As mentioned above, we can deal with them using Sqrt[1 - Cos[Indexed[\[Theta], i_]]^2] :> Sin[Indexed[\[Theta], i]]:

int[Expand@
TrigToExp[
Den1[2, 1, 1, 0]*Sin[Indexed[\[Theta], 1]]*
Sin[Indexed[\[Theta], 2]] /.
Sqrt[1 - Cos[Indexed[\[Theta], i_]]^2] :>
Sin[Indexed[\[Theta], i]]]]


We need some way of separating the exponentials, and since they're automatically combined every time you have a product, we replace them with exp[x]: /. E^(a_ + b_: 0) :> exp[a] exp[b]. Now we're nearly there:

int[Expand@
TrigToExp[
Den1[2, 1, 1, 0]*Sin[Indexed[\[Theta], 1]]*
Sin[Indexed[\[Theta], 2]] /.
Sqrt[1 - Cos[Indexed[\[Theta], i_]]^2] :>
Sin[Indexed[\[Theta], i]]] /. E^(a_ + b_: 0) :> exp[a] exp[b]]


We just need to define int for exponential terms:

int[exp[a_?(FreeQ@Indexed[_, i_]) Indexed[\[Theta], i_]], i_] =
Integrate[
Exp[a Indexed[\[Theta], i]], {Indexed[\[Theta], i], Indexed[x, i],
Pi}, {Indexed[\[Phi], i], 0, 2*Pi}]


This applies to expressions of the form $e^{a\cdot\theta_i}$, where we're using our "custom" exponential exp. Cleaning up with Simplify and ExpToTrig, we arrive at our final result (we also replace exp with Exp again):

Simplify@ExpToTrig[
int[Expand@
TrigToExp[
Den1[2, 1, 1, 0]*Sin[Indexed[\[Theta], 1]]*
Sin[Indexed[\[Theta], 2]] /.
Sqrt[1 - Cos[Indexed[\[Theta], i_]]^2] :>
Sin[Indexed[\[Theta], i]]] /.
E^(a_ + b_: 0) :> exp[a] exp[b]] /. exp -> Exp]