\( (*\ 
      This\ file\ consists\ of\ four\ blocks . \ Last\ free\ blocks\ 
        correspond\ to\ three\ parts\ of\ 
		the\ homework\ 
        problem . \ In\ the\ first\ block\ functions\ and\ parameters, \ 
      common\ to\ three\ 
		other\ blocks, \ 
      are\ defined . \ Before\ using\ three\ last\ blocks, \ 
      one\ has\ to\ execute
		the\ first\ block\ once, \ 
      to\ initilize\ those\ \(functions . \)\ *) 
	\ 
    
 (*\ n\  - \ 
        determines\ number\ of\ Matsubara\ freguencies\ which\ are\ used\ *) 
      
		
init0 := Do[
			{n = 500; 
		nn = 2\ n + 1; }]; 


    init1 := Do[

        			{init0; beta = 1/temp; 
			i1 = I; initialization}]; 


    init2 := Do[

        			{init0; temp = wmax/\((Pi\ n)\); beta = 1/temp; 
			i1 = 1; 
          initialization}]; 


    init3 := Do[

        			{init0; beta = 1/temp; 
			temp = wmax/\((Pi\ n)\); 
			
          i1 = 1; initialization}]; 	

    
 (*\ ftr[l]\  - \ 
        DFT\ corresponding\ to\ Matsubara\ frequencies
			invftr[l]\  - \ 
        inverse\ DFT
			dist[l1, l2]\  - \ norm, \ 
      which\ we\ use\ to\ check\ convergence\ 
		\ wn\  - \ 
        list\ of\ Matsubara\ frequencies
			tn\  - \ 
        list\ of\ discretized\ time\ domain
			f[x]\  - \ 
        Fermi\ function
			plott[l]\  - \ 
        plots\ functions\ in\ time\ domain
			plotf[l, x]\  - \ 
        plots\ functions\ in\ frequency\ domain, \ 
      from\  - x\ to\ x
			g0f\  - \ free\ Green' s\ function, \ 
      freq . \ representation
		\ corrf\  - 
        term\ wich\ is\ subtracted\ from\ Green' s\ function\ to\ make\ sums\ 
          covergent
			corrt\  - \ FT\ of\ corrf\ \((stepfunction)\)\ *) 
      
		

initialization := 
      Do[
		{
          t = Sqrt[1/temp]\ nn^\((\(-1\)/4)\)\ 
              Table[Exp[2\ Pi\ I\ \((n^2/2 - n\ \((s - 1)\))\)/nn], {s, 1, 
                  nn}]; 
ftr[x_] := t\ Fourier[t\ x]; 

          invftr[x_] := InverseFourier[x/t]/t; 

          dist[x_, y_] := Sum[Abs[x[\([i]\)] - y[\([i]\)]], {i, 1, nn}]/nn; 

          wn = Table[Pi\ temp\ i, {i, \(-n\), n}]; 

          tn = Table[i\ 2\ beta/nn, {i, \(-n\), n}]; 

          f[x_] := 1/\((Exp[wn[\([x + n + 1]\)]\ beta] + 1)\); 

          plott[x_] := 
            ListPlot[
              Table[{tn[\([i]\)], Chop[x[\([i]\)], 10^\((\(-8\))\)]}, {i, 1, 
                  nn}], PlotRange -> All, PlotJoined -> True]; 

          plotf[x_, wlim_] := 
            Do[nlim = Floor[wlim/\((Pi\ temp)\)]; 
              ListPlot[
                Table[{wn[\([i]\)], Chop[x[\([i]\)], 10^\((\(-8\))\)]}, {i, 
                    n - nlim + Mod[nlim, 2], n + nlim, 2}], PlotRange -> All, 
                PlotJoined -> True]; ]; 
	del =  .0000000001; 
\ 
          g0f = Table[
              1/\((i1\ wn[\([i + n + 1]\)] + I\ Sign[wn[\([i + n + 1]\)]] + 
                        If[i == 0, 1, 0])\)\ Mod[i, 2] + i1\ del, {i, \(-n\), 
                n}]; \ 

          corrf = Table[
              1/\((i1\ wn[\([i + n + 1]\)] + If[i == 0, 1, 0])\)\ Mod[i, 2] + 
                i1\ del, {i, \(-n\), n}]; 

          corrt = \(-\((Im[i1] + Re[i1] I)\)\)/2\ 
              Table[Sign[tn[\([i]\)]], {i, 1, nn}]; }]
	
    
 (*\ this\ function\ iterates\ loop\ of\ Greens\ function\ \((gf)\)\ 
        computations, \ until\ either\ covergence
			is\ reached, \ 
      or\ maximum\ number\ of\ iterations\ \((maxiter)\)\ is\ 
        met . \ Definitions\ of\ some\ 
        \(variables : 
				
          gf\  - \ Green' s\ function\ freq . \ rep . 
				gt\  - \ 
            Green' s\ function\), \ 
      time\ rep . 
				sef\  - \ selfenergy, \ 
      freq . \ rep . 
				set\  - \ selfenergy, \ 
      time\ rep . 
				g1t\  - \ gt\ from\ previous\ iteration\ *) 			
    
				
	\)\(\(computation := 
      Do[
		{gf = g0f; 
	gt = invftr[gf - corrf] + corrt; 
	iter = 0; 
          
	g1t = gt/10; 
	
          While[\((dist[gt, g1t] > eps)\) && \((\(iter++\) < maxiter)\), 

            		{
\ \ \ set = \(-i1^2\)\ u^2\ gt^3; 
\ \ \ sef = ftr[set]; 
              
\ \ \ gf = 1/\((1/g0f - sef)\); 
	\ g1t = gt; 
\ \ \ 
              gt = invftr[gf - corrf] + corrt; 
		\ 
              Print["convergence=", dist[g1t, gt]]; 
	}]; }]; \)


    \)

\(\( (*\ 
      In\ this\ part\ we\ compute\ finite\ T\ Green' s\ function\ 
          gt . \ 
				u\  - \ interaction\ U
				temp\  - \ 
        temperature\ T
				esp\  - \ 
        we\ check\ convergence\ by\ comparing\ norm\ dist[gt, gt1]\ to\ eps

          				maxiter\  - \ 
        maximum\ allowed\ number\ of\ \(iterations . \)\ *) 

u = 7; 

    temp = 1; 
eps =  .001; 
maxiter = 15; 

init1; 

computation; 


    Print["Temperature=", temp, ",", \ "U=", u]; 
plott[gt]; 

    \)\)

\(\( (*\ 
      In\ this\ part\ we\ compute\ T = 
        0\ spectral\ function . \ 
				wmax\  - \ 
          defines\ the\ range\ of\ integration\ in\ freq . \ space\ *) 


    u = 7; 
wmax = 20; 
eps =  .001; 
maxiter = 1; 

init2; 


    computation; 


Print["Zero T,", \ "U=", u]; 

    spectrf = Table[\(-Im[gf[\([i]\)]]\)\ Sign[wn[\([i]\)]], {i, 1, nn}]; 

    plotf[spectrf, 10]; 
\)\)

\(\( (*\ 
      In\ this\ part\ we\ compute\ spectral\ function\ at\ finite\ T\ , \ 
      see\ explanation\ in\ the\ solution\ *) 
			
u = 7; 
temp =  .01; 
    
wmax = 20; 
eps =  .001; 
maxiter = 1; 

init3; 

gf = g0f; 
	
    iter = 0; 
	g1f = gf/10; 
	
    While[\((dist[gf, g1f] > eps)\) && \((\(iter++\) < maxiter)\), 

      		{
	a1f = 
          Table[\(-Im[gf[\([i + n + 1]\)]\ Sign[i]]\), {i, \(-n\), n}]; 
	
        a2f = Table[a1f[\([i + n + 1]\)]\ f[i], {i, \(-n\), n}]; 
	
        a3f = Table[a1f[\([i + n + 1]\)]\ f[\(-i\)], {i, \(-n\), n}]; 
	
        c1f = Table[Re[gf[\([i + n + 1]\)]], {i, \(-n\), n}]; 
	
        a1t = invftr[a1f]; 
	a2t = invftr[a2f]; 
	a3t = invftr[a3f]; 
	
        c1t = invftr[c1f - corrf] + corrt; 
	
        imsef = 4\ u^2\ 
            ftr[a1t\ Reverse[a2t]\ a3t + Reverse[a1t]\ a2t\ a2t - 
                a1t\ Reverse[a2t]\ a2t]; 
	
        resef = 4\ u^2\ 
            ftr[c1t\ Reverse[a2t]\ a3t - Reverse[c1t]\ a2t\ a2t - 
                c1t\ Reverse[a2t]\ a2t]; 
	
        sef = Table[
            resef[\([i]\)] + I\ imsef[\([i]\)]\ Sign[\(-wn[\([i]\)]\)], {i, 
              1, nn}]; 
	g1f = gf; 
	gf = 1/\((1/g0f - sef)\); 	
	\ 
        Print["convergence=", dist[g1f, gf]]; }]; 


    Print["beta=", beta, ",", \ "U=", u]; 

    spectrf = Table[\(-Im[gf[\([i]\)]]\)\ Sign[wn[\([i]\)]], {i, 1, nn}]; 

    plotf[spectrf, 10]; 

\)\)