【托伯利兹矩阵求逆的特兰特方法】

Yasen.M

于 2022-09-16 15:53:38 发布

阅读量497

点赞数

分类专栏：常用算法（C++描述）文章标签：矩阵线性代数

本文链接：https://blog.csdn.net/qq_45449625/article/details/126891826

版权

常用算法（C++描述）专栏收录该内容

13 篇文章 0 订阅

订阅专栏

【方法说明】
设 n 阶托伯利兹矩阵为

$\boldsymbol{T}^{(n)}=\left[\begin{array}{ccccc} t_{0} & t_{1} & t_{2} & \cdots & t_{n-1} \\ \tau_{1} & t_{0} & t_{1} & \cdots & t_{n-2} \\ \tau_{2} & \tau_{1} & t_{0} & \cdots & t_{n-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \tau_{n-1} & \tau_{n-2} & \tau_{n-3} & \cdots & t_{0} \end{array}\right]$

计算程序如下：

//Toeplitz矩阵求逆.cpp
#include<cmath>
#include <iostream> 
using namespace std;
int trch (double t[],double tt[],int n,double b[])
	{ 
		int i,j,k;
		double a,s,*c,*r,*p; 
		c=new double[n]; 
		r=new double[n]; 
		p=new double[n];
		if(fabs(t[0])+1.0==1.0)
			{ 
			delete[] c;
			delete[] r;
			delete[]p; 
			cout <<"fail\n"; 
			return(0);
			}
		a=t[0];
		c[0]=tt[1]/t[0];
		r[0]=t[1]/t[0]; 
		for(k=0;k<=n-3; k++)
			{ s=0.0;
			for (j=1; j<=k+1; j++) s=s+c[k+1-j]*tt[j]; 
			s=(s-tt[k+2])/a; 
			for (i=0; i<=k; i++) p[i]=c[i]+s*r[k-i]; 
			c[k+1]=-s; 
			s=0.0;
			for (j=1; j<=k+1; j++) s=s+r[k+1-j]*t[j]; 
			s=(s-t[k+2])/a;
			for (i=0; i<=k; i++)
				{
				r[i]=r[i]+s*c[k-i]; 
				c[k-i]=p[k-i];
				}
			r[k+1]=-s; 
			a=0.0;
			for (j=1; j<=k+2; j++) a=a+t[j]*c[j-1];
			a=t[0]-a;
			if(fabs(a)+1.0==1.0)
				{ 
				delete[] c; 
				delete[] r;
				delete[] p; 
				cout <<"fail\n"; return(0);
				}
			}
		b[0]=1.0/a;
		for (i=0; i<=n-2;i++)
			{
			k=i+1;j=(i+1)*n;
			b[k]=-r[i]/a;
			b[j]=-c[i]/a;
			}
		
		for (i=0; i<=n-2; i++) 
		for (j=0; j<=n-2; j++)
			{
			k=(i+1)*n+j+1;
			b[k]=b[i*n+j]-c[i]*b[j+1]; 
			b[k]=b[k]+c[n-j-2]*b[n-i-1];
			}
		delete[] c; 
		delete[] r; 
		delete[]p; 
		return(1);
	}

int main()
{ 
	int n,i,j,k;
	double t[6]={6,5.0,4.0,3.0,2.0,1.0};
	double tt[6]={0.0,-1.0,-2.0,-3.0,-4.0,-5.0}; 
	double b[6][6],a[6][6]; 
	n=6;
	i=trch(t,tt,n,&b[0][0]); 
	if(i>0)
	{
		cout <<"B= inv(T) :\n"; 
		for(i=0;i<=5;i++)
		{
			for (j=0;j<=5;j++) 
			cout<<b[i][i] <<"\t"; 
			cout <<endl;
		}
		cout <<"A=T* B:\n"; 
		for (i=1;i<=6; i++) 
		for (j=1; j<=6; j++)
			{
			a[i-1][j-1]=0.0;
			for(k=1;k<=j-1;k++)
				a[i-1][j-1]=a[i-1][j-1]+b[i-1][k-1]*t[j-k]; 
				a[i-1][j-1]=a[i-1][j-1]+b[i-1][j-1]*t[0]; 
			for (k=j+1;k<=6;k++)
				a[i-1][j-1]=a[i-1][j-1]+b[i-1][k-1]*tt[k-j];
			}
		for (i=0; i<=5; i++)
			{
			for (j=0; j<=5;j++) cout <<a[i][j]<<"   \t"; 
			cout <<endl;
			}		
	}
	return 0;
}

计算结果: