A general implicit finite difference method for the heat equation

We study numerical methods for solving the heat equation (Subscript[u, t])=(Subscript[u, xx]) using space steps Δx and time steps Δt,and recurrences of the form

(1+2(Subscript[a, 0]))u(x,t)-(Subscript[a, 0])(u(x-Δx,t)+ u(x+Δx,t))= (1-2(∑_(j = 1)^k) (Subscript[a, j]))u(x,t-Δt)+(∑_(j = 1)^k)(Subscript[a, j]) (u(x-j Δx,t - Δt)+ u(x+j Δx,t - Δt))

for parameters (Subscript[a, 0]),(Subscript[a, 1]),...,(Subscript[a, k]). To defer consideration of boundary issues, we study the case where the space coordinate x wraps 
          around in a circle, with a period of n steps, for n even .  Our recurrence then takes the general formA u (t) = B u(t - Δt). For example, if  k=2 and  n=6  we have

k  =  2; n  = 6;
Aij[i_, j_, n_] := Switch[Mod[j - i, n], 0,(1+2(Subscript[a, 0])),1,-(Subscript[a, 0]),n - 1,-(Subscript[a, 0]),_, 0];
A[n_]  :=  Table[Aij[i, j, n] , {i,  n},  {j,  n}];
Ax  =  A[n];
Print["A = ", MatrixForm[Ax]]

and

Bij[i_, j_, k_, n_] :=
	Module[{l},l = Min[Mod[i - j, n], Mod[j - i, n]];
	If[i == j,(1-2(∑_(l = 1)^k) (Subscript[a, l])),If[l ≤ k, (Subscript[a, l]), 0]]];
B[k_, n_]  :=  Table[Bij[i, j, k, n] , {i,  n},  {j,  n}];
Bx  =  B[k, n];
Print["B = ", MatrixForm[Bx]]

Both of these matrices are recognized as convolutions .  In 
          general,the Fourier transform turns convolution into ordinary 
          multiplication .  In our discrete setting,this means that we can choose a coordinate system in which
           all convolutions can be represented as diagonal 
          matrices .  In other words, all convolution matrices can be simultaneously diagonalized
          , with respect to a coordinate system given by sin and cosine functions whose periods evenly divide our interval of n increments of Δx each .

\(\(basis[n_]\  := \ 
      Table[Cos[
          2  Floor[j/2] i\ π/n + If[j == 1, 0, Mod[j, 2]]\ π/2], 
        \ {i, n}, \ {j, n}]; 
Cx\  = \ basis[n]; \ Cxi\  = \ Inverse[Cx]; 

    Print["C = ", MatrixForm[Cx]]; 

    Print[\*""C A \!\(C\^\(-1\)\) = "", 
      MatrixForm[FullSimplify[Cxi . Ax . Cx]]]; 

    Print[\*""C B \!\(C\^\(-1\)\) = "", 
      MatrixForm[FullSimplify[Cxi . Bx . Cx]]]
\)\)

In general, we can solve for these eigenvalues by applying our recurrences to u = Cos[ϕ]using increments of  Δϕ, and then substituting  Δϕ = ((2 π m)/n) for m = 1,...,(n/2) :

(Subscript[λ, 0]) Cos[ϕ]= (1+2(Subscript[a, 0]))Cos[ϕ]-(Subscript[a, 0])(Cos[ϕ -  Δϕ]+ Cos[ϕ +  Δϕ])(Subscript[λ, 0]) Cos[ϕ]=(1+2(Subscript[a, 0]))Cos[ϕ]-2(Subscript[a, 0])Cos[ϕ]Cos[Δϕ](Subscript[λ, 0])=1+2(Subscript[a, 0])(1-Cos[Δϕ])
(Subscript[λ, 1]) Cos[ϕ]=(1-2(∑_(j = 1)^k) (Subscript[a, j]))Cos[ϕ]+(∑_(j = 1)^k)(Subscript[a, j])(Cos[ϕ -  j Δϕ]+ Cos[ϕ +  j Δϕ])
(Subscript[λ, 1]) Cos[ϕ]=(1-2(∑_(j = 1)^k) (Subscript[a, j]))Cos[ϕ]+\(\(2 ∑\)\+\(j = 1\)\%k\)(Subscript[a, j])Cos[ϕ]Cos[j Δϕ]
(Subscript[λ, 1])=1-2(∑_(j = 1)^k)(Subscript[a, j])(1-Cos[j Δϕ])
λ = ((Subscript[λ, 1])/(Subscript[λ, 0]))=((1-2(∑_(j = 1)^k)(Subscript[a, j])(1-Cos[j Δϕ]))/(1+2(Subscript[a, 0])(1-Cos[Δϕ])))

\(eigenval[k_, Δϕ_]\  := 
      \((1 - \ ∑\+\(j = 1\)\%k 2\ \(a\_j\) 
                \((1\  - \ Cos[j\ Δϕ])\))\)/
        \((1 + \ 2\ \(a\_0\) \((1\  - \ Cos[Δϕ])\))\); 

    Print[eigenval[3, Δϕ]]; 

    Table[eigenval[k, 2  π\ m/n], \ {m, n/2}]\)

On the other hand, what should these eigenvalues be? We want to study the heat equation for the initial condition u(x,0) = Cos[(ω/Δ)xx], choosing ω = ((2  π m)/( n )) so Cos[(ω/Δ)xx] completes m cycles in the interval n Δx corresponding to our circular grid. This has the general solution

\(TraditionalForm
    \`u(x, t)\  = \ 
      E\^\(\(-\((ω\/Δx)\)\^2\) t\)\ 
        Cos[ω\ Δx\ x]\)

In particular, for a time increment of Δt we have the exact solution

u(x,Δt) = (E^(-(w^2) p)) u(x, 0)  where  p = ((Δt)/(Δ(x^2)))

In other words, the desired eigenvalues of our recurrence are (E^(-(w^2) p)) for m  =  1,  ... , (n/2). We have Δϕ = ω .
For specific values of m we want choose our parameters (Subscript[a, 0]), (Subscript[a, 1]),  ... , (Subscript[a, k])so

(E^(-(w^2) p)) =  ((1-2(∑_(j = 1)^k)(Subscript[a, j])(1-Cos[jω]))/(1+2(Subscript[a, 0])(1-Cos[ω])))

or in other words we want solve a system of equations, each of the form

2  (E^(-(w^2) p))(1-Cos[ω])   (Subscript[a, 0]) + 2 (∑_(j = 1)^k)(1-Cos[jω])  (Subscript[a, j]) = 1 - (E^(-(w^2) p))    for    p = ((Δt)/(Δ(x^2))), ω = ((2  π m)/( n )).

For compatibility with other methods, in all that follows we rewrite each  (Subscript[a, j])as   (Subscript[a, j])p to make the dependence on p explicit .

idealcoefs[n_, p_, m_] :=
	Module[
		{a, b, i, j, k, x, z},
		k  =  Length[m] - 1;
		x  =  
              Table[N[E^
                    (-p 
                        (2  π 
                              m[[i]]/n)
                            ^2)],  {i, k + 1}];
		\(a\  = \ 
              Table[
				
                If[j == 1, 
					
                  x[[i]] 
                    N[2  p \((
                          1 - Cos[
                              2  π\ 
                                m[[i]] 
                                j/n])\)], 
					
                  N[2  p \((
                        1 - Cos[
                            2  π\ 
                              m[[i]] 
                              \((j - 1)\)/n])\)]], 

                				{i, k + 1}, {j, k + 1}]\);
		b  =  1 - x;
		LinearSolve[N[a], N[b]]
		]

\(desired[n_, p_, k_]\  := 
	
      Table[N[E^\((\(-p\) \((2  π\ m/n)\)^2)\)], \ {m, k}]\)

actual[n_, p_, k_, a_] :=
	Table[N[(1-2((∑_(j = 2)^(Length[a]))a[[j]] p 
                            (1 - Cos[2 π m (j - 1)/n])))/(1 + 2  a[[1]] p 
                    (1 - Cos[2  π m/n]))],
		{m, k}]

<< Graphics`MultipleListPlot`

\(graphbasic[n_, p_, k_]\  := 
	
      Module[
		{a, b, c, d, f}, 
		
        a\  = \ {desired[n, p, k], \ "true"}; 
		
        f[z_]\  := \ {actual[n, p, k, z], \ z}; 
		b\  = \ f[{0, 1}]; 

        		c\  = \ f[{1, 0}]; 
		d\  = \ f[{ .5,  .5}]; 
		
        MultipleListPlot[
			a[[1]], 
          \ b[[1]], \ 
          c[[1]], \ 
          d[[1]], 
			
          PlotJoined\  -> \ True, \ SymbolShape\  -> \ None, 
          AspectRatio\  -> 1.5, 
			
          PlotLegend\  -> 
            \ {"true", "explicit", \ "implicit", "C-N"}]; 

        	]\)

\(graphbasicerror[n_, p_, k_, t_]\  := 
	
      Module[
		{a, b, c, d, f}, 
		a\  = \ desired[n, p, k]; 
		
        f[z_]\  := \ {actual[n, p, k, z]^t - a^t, \ z}; 
		
        b\  = \ f[{0, 1}]; 
		c\  = \ f[{1, 0}]; 
		
        d\  = \ f[{ .5,  .5}]; 
		
        MultipleListPlot[
			b[[1]], 
          \ c[[1]], \ 
          d[[1]], 
			
          PlotJoined\  -> \ True, \ SymbolShape\  -> \ None, 
          AspectRatio\  -> 1.5, PlotRange\  -> \ All, 
          AxesOrigin\  -> \ {1, 0}, 
			
          PlotLegend\  -> \ {"explicit", \ "implicit", "C-N"}]; 

        	]\)

graphbasic[50,  .4, 25]

graphbasicerror[50,  .4, 10, 20]

\(\(graphsub[f_, n_, p_, m2_, m3_, m4_, k_]\  := 
	
      Module[
		{z2, z3, z4, y}, 
		z2\  = \ idealcoefs[n, p, m2]; 

        		z3\  = \ idealcoefs[n, p, m3]; 
		
        z4\  = \ idealcoefs[n, p, m4]; 

        		{f[{ .5,  .5, 0, 0}], f[z2], f[z3], f[z4]}
]	\)\)

\(graphideal[n_, p_, m2_, m3_, m4_, k_]\  := 
	
      Module[
		{a, b, c, d, e, f}, 
		
        a\  = \ {desired[n, p, k], \ "true"}; 
		
        f[z_]\  := \ {actual[n, p, k, z], \ z}; 

        		{b, c, d, e}\  = \ graphsub[f, n, p, m2, m3, m4, k]\ ; 
		
        MultipleListPlot[
			a[[1]], 
          \ b[[1]], \ 
          c[[1]], \ 
          d[[1]], \ 
          e[[1]], 
			
          PlotJoined\  -> \ True, \ SymbolShape\  -> \ None, 
          AspectRatio\  -> 1, 
			
          PlotLegend\  -> 
            \ {"true", \ "C-N", \ "2nd", \ "3rd", \ 
              "4th"}]; 
	]\)

\(grapherror[n_, p_, m2_, m3_, m4_, k_, t_]\  := 
	
      Module[
		{a, b, e, f, z2, z3}, 
		e\  = desired[n, p, k]; 
		
        f[z_]\  := \ {actual[n, p, k, z]^t - e^t, \ z}; 

        		{a, b, c, d}\  = \ graphsub[f, n, p, m2, m3, m4, k]\ ; 
		
        ListPlot[a[[1]], 
			
          PlotJoined\  -> \ True, ImageSize\  -> \ {244, 144}, 
          AspectRatio\  -> \  .5, \ PlotRange\  -> \ All, 
          AxesOrigin\  -> \ {1, 0}, \ 
          PlotLabel\  -> N[a[[2]], 4]]; 
        
		ListPlot[b[[1]], 
			
          PlotJoined\  -> \ True, \ ImageSize\  -> \ {244, 144}, 
          AspectRatio\  -> \  .5, PlotRange\  -> \ All, 
          AxesOrigin\  -> \ {1, 0}, \ 
          PlotLabel\  -> N[b[[2]], 4]]; 
        
		ListPlot[c[[1]], 
			
          PlotJoined\  -> \ True, \ ImageSize\  -> \ {244, 144}, 
          AspectRatio\  -> \  .5, PlotRange\  -> \ All, 
          AxesOrigin\  -> \ {1, 0}, \ 
          PlotLabel\  -> N[c[[2]], 4]]; 
        
		ListPlot[d[[1]], 
			
          PlotJoined\  -> \ True, ImageSize\  -> \ {244, 144}, 
          AspectRatio\  -> \  .5, \ PlotRange\  -> \ All, 
          AxesOrigin\  -> \ {1, 0}, \ 
          PlotLabel\  -> N[d[[2]], 4]]; 
        
		]\)

graphideal[50,  .5, {1, 4}, {1, 3, 5}, {1, 3, 5, 7}, 25]

grapherror[50,  .5, {1, 4}, {1, 3, 5}, {1, 3, 5, 7}, 10, 20]

graphideal[100, 2, {1, 3}, {1.5, 3.5, 25}, {1, 2.7, 4.5, 35}, 50]

grapherror[100, 2, {1, 3}, {1.5, 3.5, 25}, {1, 2.7, 4.5, 35}, 10, 30]