📄 matrix.cs
字号:
* @param eps - 计算精度
* @return bool型，求解是否成功
*/
public bool ComputeEvSymTri(double[] dblB, double[] dblC, Matrix mtxQ, int nMaxIt, double eps)
{
int i,j,k,m,it,u,v;
double d,f,h,g,p,r,e,s;
// 初值
int n = mtxQ.GetNumColumns();
dblC[n-1]=0.0;
d=0.0;
f=0.0;
// 迭代计算
for (j=0; j<=n-1; j++)
{
it=0;
h=eps*(Math.Abs(dblB[j])+Math.Abs(dblC[j]));
if (h>d)
d=h;
m=j;
while ((m<=n-1) && (Math.Abs(dblC[m])>d))
m=m+1;
if (m!=j)
{
do
{
if (it==nMaxIt)
return false;
it=it+1;
g=dblB[j];
p=(dblB[j+1]-g)/(2.0*dblC[j]);
r=Math.Sqrt(p*p+1.0);
if (p>=0.0)
dblB[j]=dblC[j]/(p+r);
else
dblB[j]=dblC[j]/(p-r);
h=g-dblB[j];
for (i=j+1; i<=n-1; i++)
dblB[i]=dblB[i]-h;
f=f+h;
p=dblB[m];
e=1.0;
s=0.0;
for (i=m-1; i>=j; i--)
{
g=e*dblC[i];
h=e*p;
if (Math.Abs(p)>=Math.Abs(dblC[i]))
{
e=dblC[i]/p;
r=Math.Sqrt(e*e+1.0);
dblC[i+1]=s*p*r;
s=e/r;
e=1.0/r;
}
else
{
e=p/dblC[i];
r=Math.Sqrt(e*e+1.0);
dblC[i+1]=s*dblC[i]*r;
s=1.0/r;
e=e/r;
}
p=e*dblB[i]-s*g;
dblB[i+1]=h+s*(e*g+s*dblB[i]);
for (k=0; k<=n-1; k++)
{
u=k*n+i+1;
v=u-1;
h=mtxQ.elements[u];
mtxQ.elements[u]=s*mtxQ.elements[v]+e*h;
mtxQ.elements[v]=e*mtxQ.elements[v]-s*h;
}
}
dblC[j]=s*p;
dblB[j]=e*p;
} while (Math.Abs(dblC[j])>d);
}
dblB[j]=dblB[j]+f;
}
for (i=0; i<=n-1; i++)
{
k=i;
p=dblB[i];
if (i+1<=n-1)
{
j=i+1;
while ((j<=n-1) && (dblB[j]<=p))
{
k=j;
p=dblB[j];
j=j+1;
}
}
if (k!=i)
{
dblB[k]=dblB[i];
dblB[i]=p;
for (j=0; j<=n-1; j++)
{
u=j*n+i;
v=j*n+k;
p=mtxQ.elements[u];
mtxQ.elements[u]=mtxQ.elements[v];
mtxQ.elements[v]=p;
}
}
}
return true;
}
/**
* 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
*/
public void MakeHberg()
{
int i = 0,j,k,u,v;
double d,t;
for (k=1; k<=numColumns-2; k++)
{
d=0.0;
for (j=k; j<=numColumns-1; j++)
{
u=j*numColumns+k-1;
t=elements[u];
if (Math.Abs(t)>Math.Abs(d))
{
d=t;
i=j;
}
}
if (d != 0.0)
{
if (i!=k)
{
for (j=k-1; j<=numColumns-1; j++)
{
u=i*numColumns+j;
v=k*numColumns+j;
t=elements[u];
elements[u]=elements[v];
elements[v]=t;
}
for (j=0; j<=numColumns-1; j++)
{
u=j*numColumns+i;
v=j*numColumns+k;
t=elements[u];
elements[u]=elements[v];
elements[v]=t;
}
}
for (i=k+1; i<=numColumns-1; i++)
{
u=i*numColumns+k-1;
t=elements[u]/d;
elements[u]=0.0;
for (j=k; j<=numColumns-1; j++)
{
v=i*numColumns+j;
elements[v]=elements[v]-t*elements[k*numColumns+j];
}
for (j=0; j<=numColumns-1; j++)
{
v=j*numColumns+k;
elements[v]=elements[v]+t*elements[j*numColumns+i];
}
}
}
}
}
/**
* 求赫申伯格矩阵全部特征值的QR方法
*
* @param dblU - 一维数组，长度为矩阵的阶数，返回时存放特征值的实部
* @param dblV - 一维数组，长度为矩阵的阶数，返回时存放特征值的虚部
* @param nMaxIt - 迭代次数
* @param eps - 计算精度
* @return bool型，求解是否成功
*/
public bool ComputeEvHBerg(double[] dblU, double[] dblV, int nMaxIt, double eps)
{
int m,it,i,j,k,l,ii,jj,kk,ll;
double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
int n = numColumns;
it=0;
m=n;
while (m!=0)
{
l=m-1;
while ((l>0) && (Math.Abs(elements[l*n+l-1]) >
eps*(Math.Abs(elements[(l-1)*n+l-1])+Math.Abs(elements[l*n+l]))))
l=l-1;
ii=(m-1)*n+m-1;
jj=(m-1)*n+m-2;
kk=(m-2)*n+m-1;
ll=(m-2)*n+m-2;
if (l==m-1)
{
dblU[m-1]=elements[(m-1)*n+m-1];
dblV[m-1]=0.0;
m=m-1;
it=0;
}
else if (l==m-2)
{
b=-(elements[ii]+elements[ll]);
c=elements[ii]*elements[ll]-elements[jj]*elements[kk];
w=b*b-4.0*c;
y=Math.Sqrt(Math.Abs(w));
if (w>0.0)
{
xy=1.0;
if (b<0.0)
xy=-1.0;
dblU[m-1]=(-b-xy*y)/2.0;
dblU[m-2]=c/dblU[m-1];
dblV[m-1]=0.0; dblV[m-2]=0.0;
}
else
{
dblU[m-1]=-b/2.0;
dblU[m-2]=dblU[m-1];
dblV[m-1]=y/2.0;
dblV[m-2]=-dblV[m-1];
}
m=m-2;
it=0;
}
else
{
if (it>=nMaxIt)
return false;
it=it+1;
for (j=l+2; j<=m-1; j++)
elements[j*n+j-2]=0.0;
for (j=l+3; j<=m-1; j++)
elements[j*n+j-3]=0.0;
for (k=l; k<=m-2; k++)
{
if (k!=l)
{
p=elements[k*n+k-1];
q=elements[(k+1)*n+k-1];
r=0.0;
if (k!=m-2)
r=elements[(k+2)*n+k-1];
}
else
{
x=elements[ii]+elements[ll];
y=elements[ll]*elements[ii]-elements[kk]*elements[jj];
ii=l*n+l;
jj=l*n+l+1;
kk=(l+1)*n+l;
ll=(l+1)*n+l+1;
p=elements[ii]*(elements[ii]-x)+elements[jj]*elements[kk]+y;
q=elements[kk]*(elements[ii]+elements[ll]-x);
r=elements[kk]*elements[(l+2)*n+l+1];
}
if ((Math.Abs(p)+Math.Abs(q)+Math.Abs(r))!=0.0)
{
xy=1.0;
if (p<0.0)
xy=-1.0;
s=xy*Math.Sqrt(p*p+q*q+r*r);
if (k!=l)
elements[k*n+k-1]=-s;
e=-q/s;
f=-r/s;
x=-p/s;
y=-x-f*r/(p+s);
g=e*r/(p+s);
z=-x-e*q/(p+s);
for (j=k; j<=m-1; j++)
{
ii=k*n+j;
jj=(k+1)*n+j;
p=x*elements[ii]+e*elements[jj];
q=e*elements[ii]+y*elements[jj];
r=f*elements[ii]+g*elements[jj];
if (k!=m-2)
{
kk=(k+2)*n+j;
p=p+f*elements[kk];
q=q+g*elements[kk];
r=r+z*elements[kk];
elements[kk]=r;
}
elements[jj]=q; elements[ii]=p;
}
j=k+3;
if (j>=m-1)
j=m-1;
for (i=l; i<=j; i++)
{
ii=i*n+k;
jj=i*n+k+1;
p=x*elements[ii]+e*elements[jj];
q=e*elements[ii]+y*elements[jj];
r=f*elements[ii]+g*elements[jj];
if (k!=m-2)
{
kk=i*n+k+2;
p=p+f*elements[kk];
q=q+g*elements[kk];
r=r+z*elements[kk];
elements[kk]=r;
}
elements[jj]=q;
elements[ii]=p;
}
}
}
}
}
return true;
}
/**
* 求实对称矩阵特征值与特征向量的雅可比法
*
* @param dblEigenValue - 一维数组，长度为矩阵的阶数，返回时存放特征值
* @param mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与数组
* dblEigenValue中第j个特征值对应的特征向量
* @param nMaxIt - 迭代次数
* @param eps - 计算精度
* @return bool型，求解是否成功
*/
public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, int nMaxIt, double eps)
{
int i,j,p = 0,q = 0,u,w,t,s,l;
double fm,cn,sn,omega,x,y,d;
if (! mtxEigenVector.Init(numColumns, numColumns))
return false;
l=1;
for (i=0; i<=numColumns-1; i++)
{
mtxEigenVector.elements[i*numColumns+i]=1.0;
for (j=0; j<=numColumns-1; j++)
if (i!=j)
mtxEigenVector.elements[i*numColumns+j]=0.0;
}
while (true)
{
fm=0.0;
for (i=1; i<=numColumns-1; i++)
{
for (j=0; j<=i-1; j++)
{
d=Math.Abs(elements[i*numColumns+j]);
if ((i!=j) && (d>fm))
{
fm=d;
p=i;
q=j;
}
}
}
if (fm<eps)
{
for (i=0; i<numColumns; ++i)
dblEigenValue[i] = GetElement(i,i);
return true;
}
if (l>nMaxIt)
return false;
l=l+1;
u=p*numColumns+q;
w=p*numColumns+p;
t=q*numColumns+p;
s=q*numColumns+q;
x=-elements[u];
y=(elements[s]-elements[w])/2.0;
omega=x/Math.Sqrt(x*x+y*y);
if (y<0.0)
omega=-omega;
sn=1.0+Math.Sqrt(1.0-omega*omega);
sn=omega/Math.Sqrt(2.0*sn);
cn=Math.Sqrt(1.0-sn*sn);
fm=elements[w];
elements[w]=fm*cn*cn+elements[s]*sn*sn+elements[u]*omega;
elements[s]=fm*sn*sn+elements[s]*cn*cn-elements[u]*omega;
elements[u]=0.0;
elements[t]=0.0;
for (j=0; j<=numColumns-1; j++)
{
if ((j!=p) && (j!=q))
{
u=p*numColumns+j; w=q*numColumns+j;
fm=elements[u];
elements[u]=fm*cn+elements[w]*sn;
elements[w]=-fm*sn+elements[w]*cn;
}
}
for (i=0; i<=numColumns-1; i++)
{
if ((i!=p) && (i!=q))
{
u=i*numColumns+p;
w=i*numColumns+q;
fm=elements[u];
elements[u]=fm*cn+elements[w]*sn;
elements[w]=-fm*sn+elements[w]*cn;
}
}
for (i=0; i<=numColumns-1; i++)
{
u=i*numColumns+p;
w=i*numColumns+q;
fm=mtxEigenVector.elements[u];
mtxEigenVector.elements[u]=fm*cn+mtxEigenVector.elements[w]*sn;
mtxEigenVector.elements[w]=-fm*sn+mtxEigenVector.elements[w]*cn;
}
}
}
/**
* 求实对称矩阵特征值与特征向量的雅可比过关法
*
* @param dblEigenValue - 一维数组，长度为矩阵的阶数，返回时存放特征值
* @param mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与数组
* dblEigenValue中第j个特征值对应的特征向量
* @param eps - 计算精度
* @return bool型，求解是否成功
*/
public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, double eps)
{
int i,j,p,q,u,w,t,s;
double ff,fm,cn,sn,omega,x,y,d;
if (! mtxEigenVector.Init(numColumns, numColumns))
return false;
for (i=0; i<=numColumns-1; i++)
{
mtxEigenVector.elements[i*numColumns+i]=1.0;
for (j=0; j<=numColumns-1; j++)
if (i!=j)
mtxEigenVector.elements[i*numColumns+j]=0.0;
}
ff=0.0;
for (i=1; i<=numColumns-1; i++)
{
for (j=0; j<=i-1; j++)
{
d=elements[i*numColumns+j];
ff=ff+d*d;
}
}
ff=Math.Sqrt(2.0*ff);
ff=ff/(1.0*numColumns);
bool nextLoop = false;
while (true)
{
for (i=1; i<=numColumns-1; i++)
{
for (j=0; j<=i-1; j++)
{
d=Math.Abs(elements[i*numColumns+j]);
if (d>ff)
{
p=i;
q=j;
u=p*numColumns+q;
w=p*numColumns+p;
t=q*numColumns+p;
s=q*numColumns+q;
x=-elements[u];
y=(elements[s]-elements[w])/2.0;
omega=x/Math.Sqrt(x*x+y*y);
if (y<0.0)
omega=-omega;
sn=1.0+Math.Sqrt(1.0-omega*omega);
sn=omega/Math.Sqrt(2.0*sn);
cn=Math.Sqrt(1.0-sn*sn);
fm=elements[w];
elements[w]=fm*cn*cn+elements[s]*sn*sn+elements[u]*omega;
elements[s]=fm*sn*sn+elements[s]*cn*cn-elements[u]*omega;
elements[u]=0.0; elements[t]=0.0;
for (j=0; j<=numColumns-1; j++)
{
if ((j!=p)&&(j!=q))
{
u=p*numColumns+j;
w=q*numColumns+j;
fm=elements[u];
elements[u]=fm*cn+elements[w]*sn;
elements[w]=-fm*sn+elements[w]*cn;
}
}
for (i=0; i<=numColumns-1; i++)
{
if ((i!=p)&&(i!=q))
{
u=i*numColumns+p;
w=i*numColumns+q;
fm=elements[u];
elements[u]=fm*cn+elements[w]*sn;
elements[w]=-fm*sn+elements[w]*cn;
}
}
for (i=0; i<=numColumns-1; i++)
{
u=i*numColumns+p;
w=i*numColumns+q;
fm=mtxEigenVector.elements[u];
mtxEigenVector.elements[u]=fm*cn+mtxEigenVector.elements[w]*sn;
mtxEigenVector.elements[w]=-fm*sn+mtxEigenVector.elements[w]*cn;
}
nextLoop = true;
break;
}
}
if (nextLoop)
break;
}
if (nextLoop)
{
nextLoop = false;
continue;
}
nextLoop = false;
// 如果达到精度要求，退出循环，返回结果
if (ff<eps)
{
for (i=0; i<numColumns; ++i)
dblEigenValue[i] = GetElement(i,i);
return true;
}
ff=ff/(1.0*numColumns);
}
}
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -