0
$\begingroup$

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?

$\endgroup$
  • 3
    $\begingroup$ 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 *) $\endgroup$ – J. M. is away Aug 24 '17 at 7:22
  • 1
    $\begingroup$ Simply putting an Evaluate around the function in the plot should help a lot. $\endgroup$ – mmeent Aug 24 '17 at 7:27
  • $\begingroup$ 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). $\endgroup$ – aardvark2012 Aug 24 '17 at 10:46
0
$\begingroup$

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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