ToeplitzSolve Routines

Summary

Solve a toepltiz system of linear equations.

Declaration

Procedure ToeplitzSolve(FirstRow, B, X: TVec);

Description

Solves toeplitz system. Solves toeplitz system, defined as:


     [  R(0)     R(1)* ...   R(N-1)* ]       [ A[1]   ]  = [  -B[1]  ]
     [  R(1)     R(0)  ...   R(N-2)* ]       [ A[2]   ]  = [  -B[2]  ]
     [   ....        .         .     ]   x   [   .    ]  = [    .    ]
     [  R(N-2)   R(N-3) ...   R(1)*  ]       [ A[N-1] ]  = [ -B[N-1] ]
     [  R(N-1)   R(N-2) ...   R(0)   ]       [ A[N]   ]  = [ -B[N]   ]

Declaration

Procedure ToeplitzSolve(FirstRow, FirstCol, B, X: TVec);

Description

The procedures solves Toeplitz system of linear equations defined as:


     [  R(0)     R(1) ...    R(N-1)  ]       [ X(0)  ]  = [  -B(0)  ]
     [  C(0)     R(0)  ...   R(N-2)  ]       [ X(1)  ]  = [  -B(1)  ]
     [   ....        .         .     ] x     [   .   ]  = [    .    ]
     [  C(N-1) C(N-2) ...    R(1)    ]       [ X(N-2)]  = [ -B(N-2) ]
     [  C(N-2) C(N-1) ...    R(0)    ]       [ X(N-1)]  = [ -B(N-1) ]

The computational complexity of the algorithm is O(n2) as oposed to at least O(n3) for solving a general matrix. The solution is stored in vector X. B, FirstRow and FirstCol can be real or complex. R stands for First Row and C stands for First column. FirstCol.Length must be equal to FirstRow.Length-1.

See Also
Levinson

Example 1

    uses MtxExpr, Math387, MtxVec, MtxVecTee, MtxVecEdit, Toeplitz;

    procedure TForm1.Button1Click(Sender: TObject);
    var FirstRow,FirstCol,X1,X2,B,R: Vector;
        A: Matrix;
    begin
         FirstRow.Size(4,False);
         FirstCol.Size(3,False);
         B.Size(4,false);
         FirstRow[0] := 1;
         FirstRow[1] := 2;
         FirstRow[2] := 3;
         FirstRow[3] := 4;
         FirstCol[0] := 2;
         FirstCol[1] := 3;
         FirstCol[2] := 4;
         A.Toeplitz(FirstRow,FirstCol);
         B[0] := -2;
         B[1] := -3;
         B[2] := -4;
         B[3] := -5;
         R.Copy(FirstRow);
         R.Resize(5);
         R[4] := 5;
         ToeplitzSolve(FirstRow,FirstCol,B,X1);
         A.LUSolve(B,X2,mtGeneral);
         if not X2.Equal(X1,EPS,cmpAbsolute) then Eraise('Not equal!');
         ViewValues(X1,'Real X1',true);
         Levinson(R,X2);
         ViewValues(X2,'Real X2',true);

         //       Complex
         FirstRow.Size(4,True);
         FirstCol.Size(3,True);
         B.Size(4,True);

         FirstRow.CValues[0] := Cplx(1,0);
         FirstRow.CValues[1] := Cplx(2,0);
         FirstRow.CValues[2] := Cplx(3,-1);
         FirstRow.CValues[3] := Cplx(4,0);
         FirstCol.CValues[0] := Cplx(2,0);
         FirstCol.CValues[1] := Cplx(3,1);
         FirstCol.CValues[2] := Cplx(4,0);
         A.Toeplitz(FirstRow,FirstCol);
         B.CValues[0] := Cplx(-2,0);
         B.CValues[1] := Cplx(-3,-1);
         B.CValues[2] := Cplx(-4,0);
         B.CValues[3] := Cplx(-5,0);
         R.Copy(FirstRow);
         R.Resize(5);
         R.CValues[4] := Cplx(5,0);
         ToeplitzSolve(FirstRow,FirstCol,B,X1);
         A.LUSolve(B,X2,mtGeneral);
         if not X2.Equal(X1,EPS*50,cmpAbsolute) then Eraise('Not equal!');
         ViewValues(X1,'Complex X1',true);
         Levinson(R,X2);
         ViewValues(X2,'Complex X2',true);
    end;

    #include "MtxVecCPP.h" //MtxVecCPP.cpp must be included in the project
    #include "MtxVecEdit.hpp"
    #include "MtxVecTee.hpp"
    #include "Toeplitz.hpp"

    void __fastcall TForm1::BitBtn1Click(TObject *Sender)
    {
        Vector FirstRow,FirstCol,X1,X2,B,R;
        Matrix	A;

        FirstRow->Size(4,false);
        FirstCol->Size(3,false);
        B->Size(4,false);
        FirstRow[0] = 1;
        FirstRow[1] = 2;
        FirstRow[2] = 3;
        FirstRow[3] = 4;
        FirstCol[0] = 2;
        FirstCol[1] = 3;
        FirstCol[2] = 4;
        A->Toeplitz(FirstRow,FirstCol);
        B[0] = -2;
        B[1] = -3;
        B[2] = -4;
        B[3] = -5;
        R->Copy(FirstRow);
        R->Resize(5);
        R[4] = 5;
        ToeplitzSolve(FirstRow,FirstCol,B,X1);
        A->LUSolve(B,X2,mtGeneral);
        if ((X2->Equal(X1,10*EPS,cmpAbsolute)) == false)
        {
          throw "Not equal";
        }
        ViewValues(X1,"Real X1",true);
        Levinson(R,X2);
        ViewValues(X2,"Real X2",true);
        //       Complex

        FirstRow->Size(4,true);
        FirstCol->Size(3,true);
        B->Size(4,True);

        FirstRow->CValues[0] = Cplx(1,0);
        FirstRow->CValues[1] = Cplx(2,0);
        FirstRow->CValues[2] = Cplx(3,-1);
        FirstRow->CValues[3] = Cplx(4,0);
        FirstCol->CValues[0] = Cplx(2,0);
        FirstCol->CValues[1] = Cplx(3,1);
        FirstCol->CValues[2] = Cplx(4,0);
        A->Toeplitz(FirstRow,FirstCol);
        B->CValues[0] = Cplx(-2,0);
        B->CValues[1] = Cplx(-3,-1);
        B->CValues[2] = Cplx(-4,0);
        B->CValues[3] = Cplx(-5,0);
        R->Copy(FirstRow);
        R->Resize(5);
        R->CValues[4] = Cplx(5,0);
        ToeplitzSolve(FirstRow,FirstCol,B,X1);
        A->LUSolve(B,X2,mtGeneral);
        if (!X2->Equal(X1,EPS*50,cmpAbsolute))
        {
          throw "Not equal!";
        }
        ViewValues(X1,"Complex X1",true);
        Levinson(R,X2);
        ViewValues(X2,"Complex X2",true);
    }

using Dew.Math;
using Dew.Math.Units;
using Dew.Math.Editors;

namespace Dew.Examples
{
	private void Example()
	{

    Vector FirstRow = new Vector(0);
		Vector FirstCol = new Vector(0);
		Vector X1 = new Vector(0);
		Vector X2 = new Vector(0);
		Vector B = new Vector(0);
		Vector R = new Vector(0;
		Matrix A = new Matrix(0,0);

    FirstRow.SetIt(false,new double[] {1,2,3,4});
    FirstCol.SetIt(false,new double[] {2,3,4});
    A.Toeplitz(FirstRow,FirstCol);
		B.SetIt(false, new double[] {-2,-3,-4,-5});
    R.Copy(FirstRow);
		R.Resize(5);
    R.Values[4] = 5;
	  Toeplitz.ToeplitzSolve(FirstRow,FirstCol,B,X1);
	  A.LUSolve(B,X2,MtxVec.mtGeneral);
    MtxVecEdit.ViewValues(X1,"Real X1",true);
   	Toeplitz.Levinson(R,X2);
   	MtxVecEdit.ViewValues(X2,"Real X2",true);
	}
}