Example 2

Let's find the eigenfunction expansion of the Heaviside function on -1<x<1.

In[23]:=

  Clear[f,c0,c,Sn,fplot]

In[24]:=

  f[x_]:=UnitStep[x]

In[25]:=

  fplot=Plot[f[x],{x,-1,1}];

Recall that

This converges to f(x) where f(x) is continuous and converges to
1/2 [f(x+)+f(x-)] where f(x) has a jump discontinuity.

Remember, with Legendre polynomials p(x)=1, the eigenfunctions are the Legendre polynomials, and the limits of integration are -1 and 1.

Recall from what we saw above that the denominator in the expression for cn is 2/(2n+1).

In[26]:=

  c0=1/2 Integrate[f[x]LegendreP[0,x],{x,-1,1}];
  Simplify[c0]
  c=Table[(2k+1)/2 Integrate[f[x]LegendreP[k,x],{x,-1,1}],{k,1,6}];
  Simplify[c]

Out[26]=

  1
  -
  2

Out[27]=

   3       7       11
  {-, 0, -(--), 0, --, 0}
   4       16      32

In[28]:=

  Sn=Table[c0 LegendreP[0,x]+Sum[c[[k]]LegendreP[k,x],{k,1,n}],{n,1,6}];

Let's see how we did.

In[29]:=

  sn1=Plot[Sn[[1]],{x,-3,3},DisplayFunction->Identity];
  Show[fplot,sn1,DisplayFunction->$DisplayFunction,
  PlotLabel->FontForm["c0P0+c1P1",{"Helvetica-Bold",12}]];

In[30]:=

  sn3=Plot[Sn[[3]],{x,-3,3},DisplayFunction->Identity];
  Show[fplot,sn3,PlotLabel->
  FontForm["c0P0+c1P1+c2P2+c3P3",{"Helvetica-Bold",12}]];

In[31]:=

  sn6=Plot[Sn[[6]],{x,-3,3},DisplayFunction->Identity];
  Show[fplot,sn6, PlotLabel->
  FontForm["Sum(n=0 to 6) cnPn",{"Helvetica-Bold",12}]];

Let's go a little further this time.

In[32]:=

  Clear[c0,c,Sn,sn]
  c0=1/2 Integrate[f[x]LegendreP[0,x],{x,-1,1}];
  Simplify[c0]
  c=Table[(2k+1)/2 Integrate[f[x]LegendreP[k,x],{x,-1,1}],{k,1,10}];
  Simplify[c]

Out[32]=

  1
  -
  2

Out[33]=

   3       7       11       75       133
  {-, 0, -(--), 0, --, 0, -(---), 0, ---, 0}
   4       16      32       256      512

In[34]:=

  Sn=Table[c0 LegendreP[0,x]+Sum[c[[k]]LegendreP[k,x],{k,1,n}],{n,1,10}];

In[35]:=

  sn8=Plot[Sn[[8]],{x,-1.5,1.5},DisplayFunction->Identity];
  Show[fplot,sn8, PlotLabel->
  FontForm["Sum(n=0 to 8) cnPn",{"Helvetica-Bold",12}]];

In[36]:=

  sn10=Plot[Sn[[10]],{x,-1.5,1.5},DisplayFunction->Identity];
  Show[fplot,sn10, PlotLabel->
  FontForm["Sum(n=0 to 10) cnPn",{"Helvetica-Bold",12}]];

You can modify this to go with as many terms as you want (or as many terms as you are willing to wait for).

Up to Legendre Polynomials