声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2908|回复: 8

[控制理论] [求助]求教大家,平方根辨识问题?

[复制链接]
发表于 2006-3-22 23:44 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
<P>近期遇到一个问题,在一篇论文上看到<br>A*X=0,其中A为已知的数据阵,X为向量.论文上说求X用平方根法辨识,<br>但不知是怎么求得?<br>在此请教各位大侠</P>
[此贴子已经被cdwxg于2006-5-28 0:32:08编辑过]

回复
分享到:

使用道具 举报

 楼主| 发表于 2006-3-23 21:07 | 显示全部楼层
<FONT color=#f73809 size=4>非常感谢yejet的解答.<br>万分感谢!!!!!!!<br>以后有问题再向你请教.</FONT>[em01][em01]
[此贴子已经被作者于2006-3-23 22:32:36编辑过]

 楼主| 发表于 2006-3-23 22:56 | 显示全部楼层
刚看完yejet提供的资料,想问这个资料摘自何处?<BR>可能是我现在一下子没看明白,所以我还想向你请教一下,如何辨识的问题,即A*X=0(其中A为n*n已知矩阵,X为一个列向量),如何求X.
发表于 2006-4-7 19:54 | 显示全部楼层

回复:(liljx_2008)刚看完yejet提供的资料,想问这个...

<DIV class=quote><B>以下是引用<I>liljx_2008</I>在2006-3-23 22:56:16的发言:</B><BR>刚看完yejet提供的资料,想问这个资料摘自何处?<BR>可能是我现在一下子没看明白,所以我还想向你请教一下,如何辨识的问题,即A*X=0(其中A为n*n已知矩阵,X为一个列向量),如何求X.</DIV>
<br>??????上面不是说得很清楚吗?
发表于 2006-4-14 16:04 | 显示全部楼层
是呀,文章中的变量说明也不全面
发表于 2006-4-15 11:13 | 显示全部楼层

回复:(yejet)回复:(liljx_2008)求教大家,平方根...

<FONT face=Verdana color=#da2549><B>yejet, 您有相应的程序吗?或者推荐一个相似的程序也行,我急用。先谢了!</B></FONT>
[此贴子已经被作者于2006-4-15 11:17:03编辑过]

发表于 2006-4-15 19:54 | 显示全部楼层

回复:(liljx_2008)求教大家,平方根辨识问题?

线性方程组的求解---平方根法及改进平方根法<BR><BR>
<P>平方根法主要用于求解对称正定矩阵方程:</P>
<P>首先要提到的是有个关于正定矩阵的定理,说的是若A为n阶地称正定矩阵,则存在一个实的非奇异下三角矩阵L,使A=LL’(L’表示L的对称矩阵)</P>
<P>根据这个前提,在结合前面的LU分解算法,便有了这里的平方根算法:</P>
<P>
<br>
<p>
<P>平方根方法:</P>
<P>/*
<p>
<p>
<P>squareRootMethod, coded by EmilMathew 05/9/10, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.
<p>
<p>
<P>*/
<p>
<p>
<P>void squareRootMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,int size)
<p>
<p>
<P>{
<p>
<p>
<P>/*Maths Reason:
<p>
<p>
<P>L*U=A*x=b
<p>
<p>
<P>When you meet a duiCheng and zheng ding matrix:
<p>
<p>
<P>it could be pation like this:
<p>
<p>
<P>l11 l11 l21 ... ln1
<p>
<p>
<P>l21 l22 * l22 ... ln2
<p>
<p>
<P>.... ... ...
<p>
<p>
<P>ln1 ln2 ... lnn lnn
<p>
<p>
<P>
<p>
<p>
<P>i=j:
<p>
<p>
<P>lij=sqrt(aij-sigma_k1...j-1(ljk^2))
<p>
<p>
<P>
<p>
<p>
<P>i&gt;j:
<p>
<p>
<P>lij=(aij-sigma_k1...j-1(lik*ljk))/ljj
<p>
<p>
<P>
<p>
<p>
<P>the steps below is very easy :
<p>
<p>
<P>L*y=b;
<p>
<p>
<P>U*x=y;
<p>
<p>
<P>
<p>
<p>
<P>Enjoy!:)
<p>
<p>
<P>by EmilMatthew
<p>
<p>
<P>05/9/10.
<p>
<p>
<P>*/
<p>
<p>
<P>Type** l_Matrix,* yAnsList;
<p>
<p>
<P>Type tmpData;
<p>
<p>
<P>int i,j;
<p>
<p>
<P>
<p>
<p>
<P>/*pointer data assertion*/
<p>
<p>
<P>assertF(inMatrixArr!=NULL,"in squareRootMethod,matrixArr is NULL\n");
<p>
<p>
<P>assertF(bList!=NULL,"in squareRootMethod,bList is NULL\n");
<p>
<p>
<P>assertF(xAnsList!=NULL,"in squareRootMethod,xAnsList is NULL\n");
<p>
<p>
<P>
<p>
<p>
<P>/*correct pass in matrix assertion*/
<p>
<p>
<P>assertF(duiChengMatrixCheck(inMatrixArr,size),"in squareRootMethod,the pass in matrix is not dui cheng\n");
<p>
<p>
<P>
<p>
<p>
<P>/*Mem Apply*/
<p>
<p>
<P>listArrMemApply(&amp;yAnsList,size);
<p>
<p>
<P>twoDArrMemApply(&amp;l_Matrix,size,size);
<p>
<p>
<P>
<p>
<p>
<P>assertF(l_Matrix!=NULL,"in squareRootMethod,l_Matrix is null\n");
<p>
<p>
<P>assertF(yAnsList!=NULL,"in squareRootMethod,yAnsList is null\n");
<p>
<p>
<P>
<p>
<p>
<P>/*Core Program*/
<p>
<p>
<P>for(i=0;i
<p>
<P>for(j=0;j&lt;=i;j++)
<p>
<p>
<P>{
<p>
<p>
<P>if(i==j)
<p>
<p>
<P>{
<p>
<p>
<P>tmpData=sumSomeRowPower(l_Matrix,j,0,j-1,2);
<p>
<p>
<P>// printf("tmpData:%f\n",tmpData);
<p>
<p>
<P>l_Matrix[j]=(float)sqrt(inMatrixArr-tmpData);
<p>
<p>
<P>}
<p>
<p>
<P>else
<p>
<p>
<P>{
<p>
<p>
<P>l_Matrix[j]=(inMatrixArr[j]-sumTwoRowBy(l_Matrix,i,j,0,j-1))/l_Matrix[j][j];
<p>
<p>
<P>}
<p>
<p>
<P>}
<p>
<p>
<P>for(i=0;i
<p>
<P>yAnsList=(bList-sumArr_JKByList_K(l_Matrix,yAnsList,i,0,i-1))/l_Matrix;
<p>
<p>
<P>
<p>
<p>
<P>for(i=size-1;i&gt;=0;i--)
<p>
<p>
<P>xAnsList=(yAnsList-sumArr_KJByList_K(l_Matrix,xAnsList,i,i+1,size-1))/l_Matrix;
<p>
<p>
<P>
<p>
<p>
<P>twoDArrMemFree(&amp;l_Matrix,size);
<p>
<p>
<P>free(yAnsList);
<p>
<p>
<P>}
<p>
<p>
<P>
<p>
<p>
<P>平方根算法的计算量约为普通三角分解算法的一半,但由于这里要用到开平方,效率不是很高,所以,便有了为效率而存在的改进版平方根算法:
<p>
<p>
<P>改进的平方根算法:
<p>
<p>
<P>/*
<p>
<p>
<P>enhancedSquareRootMethod, coded by EmilMathew 05/9/10, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.
<p>
<p>
<P>*/
<p>
<p>
<P>
<p>
<p>
<P>如下:
<p>
<p>
<P>void enhancedSquareRootMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,int size)
<p>
<p>
<P>{
<p>
<p>
<P>Type** l_Matrix,** t_Matrix;
<p>
<p>
<P>Type* yAnsList,* dList;
<p>
<p>
<P>
<p>
<p>
<P>int i,j;
<p>
<p>
<P>
<p>
<p>
<P>/*pointer data assertion*/
<p>
<p>
<P>assertF(inMatrixArr!=NULL,"in enhancedSquareRootMethod,matrixArr is NULL\n");
<p>
<p>
<P>assertF(bList!=NULL,"in enhancedSquareRootMethod,bList is NULL\n");
<p>
<p>
<P>assertF(xAnsList!=NULL,"in enhancedSquareRootMethod,xAnsList is NULL\n");
<p>
<p>
<P>
<p>
<p>
<P>/*correct pass in matrix assertion*/
<p>
<p>
<P>assertF(duiChengMatrixCheck(inMatrixArr,size),"in enhancedSquareRootMethod,the pass in matrix is not dui cheng\n");
<p>
<p>
<P>
<p>
<p>
<P>/*Mem Apply*/
<p>
<p>
<P>listArrMemApply(&amp;yAnsList,size);
<p>
<p>
<P>listArrMemApply(&amp;dList,size);
<p>
<p>
<P>twoDArrMemApply(&amp;l_Matrix,size,size);
<p>
<p>
<P>twoDArrMemApply(&amp;t_Matrix,size,size);
<p>
<p>
<P>
<p>
<p>
<P>assertF(t_Matrix!=NULL,"in enhancedSquareRootMethod,t_Matrix is null\n");
<p>
<p>
<P>assertF(l_Matrix!=NULL,"in enhancedSquareRootMethod,l_Matrix is null\n");
<p>
<p>
<P>assertF(yAnsList!=NULL,"in enhancedSquareRootMethod,yAnsList is null\n");
<p>
<p>
<P>
<p>
<p>
<P>for(i=0;i
<p>
<P>l_Matrix=1;
<p>
<p>
<P>
<p>
<p>
<P>for(j=0;j
<p>
<P>{
<p>
<p>
<P>dList[j]=inMatrixArr[j][j]-sumArr1_IKByArr2_JK(t_Matrix,l_Matrix,j,j,0,j-1);
<p>
<p>
<P>for(i=j+1;i
<p>
<P>{
<p>
<p>
<P>t_Matrix[j]=inMatrixArr[j]-sumArr1_IKByArr2_JK(inMatrixArr,l_Matrix,i,j,0,j-1);
<p>
<p>
<P>l_Matrix[j]=t_Matrix[j]/dList[j];
<p>
<p>
<P>}
<p>
<p>
<P>}
<p>
<p>
<P>
<p>
<p>
<P>
<p>
<p>
<P>for(i=0;i
<p>
<P>yAnsList=bList-sumArr_JKByList_K(l_Matrix,yAnsList,i,0,i-1);
<p>
<p>
<P>
<p>
<p>
<P>for(i=size-1;i&gt;=0;i--)
<p>
<p>
<P>xAnsList=yAnsList/dList-sumArr_KJByList_K(l_Matrix,xAnsList,i,i+1,size-1);
<p>
<p>
<P>
<p>
<p>
<P>/*mem free*/
<p>
<p>
<P>twoDArrMemFree(&amp;t_Matrix,size);
<p>
<p>
<P>twoDArrMemFree(&amp;l_Matrix,size);
<p>
<p>
<P>free(yAnsList);
<p>
<p>
<P>free(dList);
<p>
<p>
<P>}
<p>
<p>
<P>
<p>
<p>
<P>测试程序:
<p>
<p>
<P>/*Square Method Algorithm test program*/
<p>
<p>
<P>#include "Global.h"
<p>
<p>
<P>#include "Ulti.h"
<p>
<p>
<P>#include "MyAssert.h"
<p>
<p>
<P>#include "Matrix.h"
<p>
<p>
<P>#include
<p>
<p>
<P>#include
<p>
<p>
<P>#include
<p>
<p>
<P>#include
<p>
<p>
<P>#include
<p>
<p>
<P>
<p>
<p>
<P>
<p>
<p>
<P>char *inFileName="inputData.txt";
<p>
<p>
<P>/*
<p>
<p>
<P>input data specification
<p>
<p>
<P>len,
<p>
<p>
<P>a00,a01,...,a0n-1,b0;
<p>
<p>
<P>.....
<p>
<p>
<P>an-10,an-11,...,an-1n-1,bn-1;
<p>
<p>
<P>*/
<p>
<p>
<P>
<p>
<p>
<P>
<p>
<p>
<P>char *outFileName="outputData.txt";
<p>
<p>
<P>#define DEBUG 1
<p>
<p>
<P>
<p>
<p>
<P>void main(int argc,char* argv[])
<p>
<p>
<P>{
<p>
<p>
<P>FILE *inputFile;/*input file*/
<p>
<p>
<P>FILE *outputFile;/*output file*/
<p>
<p>
<P>
<p>
<p>
<P>double startTime,endTime,tweenTime;/*time callopsed info*/
<p>
<p>
<P>
<p>
<p>
<P>/*The read in data*/
<p>
<p>
<P>int len,methodIndex;
<p>
<p>
<P>Type** matrixArr;
<p>
<p>
<P>Type* bList,* xAnsList;
<p>
<p>
<P>
<p>
<p>
<P>
<p>
<p>
<P>int i,j;/*iterator index*/
<p>
<p>
<P>
<p>
<p>
<P>/*input file open*/
<p>
<p>
<P>if(argc&gt;1)strcpy(inFileName,argv[1]);
<p>
<p>
<P>assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");
<p>
<p>
<P>printf("input file open success\n");
<p>
<p>
<P>
<p>
<p>
<P>/*outpout file open*/
<p>
<p>
<P>if(argc&gt;2)strcpy(outFileName,argv[2]);
<p>
<p>
<P>assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");
<p>
<p>
<P>printf("output file open success\n");
<p>
<p>
<P>
<p>
<p>
<P>fscanf(inputFile,"size=%d;\r\n",&amp;len);
<p>
<p>
<P>fscanf(inputFile,"method=%d;\r\n",&amp;methodIndex);
<p>
<p>
<P>
<p>
<p>
<P>/*Memory apply*/
<p>
<p>
<P>matrixArr=(Type**)malloc(sizeof(Type*)*len);
<p>
<p>
<P>for(i=0;i
<p>
<P>matrixArr=(Type*)malloc(sizeof(Type)*len);
<p>
<p>
<P>
<p>
<p>
<P>bList=(Type*)malloc(sizeof(Type)*len);
<p>
<p>
<P>xAnsList=(Type*)malloc(sizeof(Type)*len);
<p>
<p>
<P>
<p>
<p>
<P>/*Read info data*/
<p>
<p>
<P>for(i=0;i
<p>
<P>{
<p>
<p>
<P>for(j=0;j
<p>
<P>fscanf(inputFile,"%f,",&amp;matrixArr[j]);
<p>
<p>
<P>fscanf(inputFile,"%f;",&amp;bList);
<p>
<p>
<P>}
<p>
<p>
<P>
<p>
<p>
<P>
<p>
<p>
<P>/*Check the input data*/
<p>
<p>
<P>showArrListFloat(bList,0,len);
<p>
<p>
<P>show2DArrFloat(matrixArr,len,len);
<p>
<p>
<P>
<p>
<p>
<P>#if DEBUG
<p>
<p>
<P>printf("\n*******start of test program******\n");
<p>
<p>
<P>printf("now is runnig,please wait...\n");
<p>
<p>
<P>startTime=(double)clock()/(double)CLOCKS_PER_SEC;
<p>
<p>
<P>/******************Core program code*************/
<p>
<p>
<P>switch(methodIndex)
<p>
<p>
<P>{
<p>
<p>
<P>case 1:
<p>
<p>
<P>enhancedSquareRootMethod(matrixArr,bList,xAnsList,len);
<p>
<p>
<P>printf("after the enhancedSquareRootMethod:the ans x rows is:\n");
<p>
<p>
<P>fprintf(outputFile,"after the enhancedSquareRootMethod:the ans x rows is:(from x0 to xn-1)\r\n");
<p>
<p>
<P>break;
<p>
<p>
<P>
<p>
<p>
<P>case 2:
<p>
<p>
<P>
<p>
<p>
<P>squareRootMethod(matrixArr,bList,xAnsList,len);
<p>
<p>
<P>printf("after the SquartRootPationMethod:the ans x rows is:\n");
<p>
<p>
<P>fprintf(outputFile,"after the SquartRootPationMethod:the ans x rows is:(from x0 to xn-1)\r\n");
<p>
<p>
<P>break;
<p>
<p>
<P>
<p>
<p>
<P>default:
<p>
<p>
<P>printf("input method index error\n");
<p>
<p>
<P>break;
<p>
<p>
<P>}
<p>
<p>
<P>showArrListFloat(xAnsList,0,len);
<p>
<p>
<P>outputListArrFloat(xAnsList,0,len,outputFile);
<p>
<p>
<P>/******************End of Core program**********/
<p>
<p>
<P>endTime=(double)clock()/(double)CLOCKS_PER_SEC;
<p>
<p>
<P>tweenTime=endTime-startTime;/*Get the time collapsed*/
<p>
<p>
<P>/*Time collapsed output*/
<p>
<p>
<P>printf("the collapsed time in this algorithm implement is:%f\n",tweenTime);
<p>
<p>
<P>fprintf(outputFile,"the collapsed time in this algorithm implement is:%f\r\n",tweenTime);
<p>
<p>
<P>printf("\n*******end of test program******\n");
<p>
<p>
<P>#endif
<p>
<p>
<P>for(i=0;i
<p>
<P>free(matrixArr);
<p>
<p>
<P>free(matrixArr);
<p>
<p>
<P>
<p>
<p>
<P>free(xAnsList);
<p>
<p>
<P>free(bList);
<p>
<p>
<P>
<p>
<p>
<P>printf("program end successfully,\n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer.\n");
<p>
<p>
<P>getchar();/*Screen Delay Control*/
<p>
<p>
<P>return;
<p>
<p>
<P>}
<p>
<p>
<P>
<p>
<p>
<P>测试结果:
<p>
<p>
<P>平方根法:
<p>
<p>
<P>test1:
<p>
<p>
<P>输入:
<p>
<p>
<P>size=3;
<p>
<p>
<P>5,-4,1,1;
<p>
<p>
<P>-4,6,-4,2;
<p>
<p>
<P>1,-4,6,-3;
<p>
<p>
<P>输出:
<p>
<p>
<P>after the SquartRootPationMethod:the ans x rows is:(from x0 to xn-1)
<p>
<p>
<P>1.00000 1.00000 0.00000
<p>
<p>
<P>
<p>
<p>
<P>改进平方根法:
<p>
<p>
<P>test1:
<p>
<p>
<P>输入:
<p>
<p>
<P>size=3;
<p>
<p>
<P>method=1;
<p>
<p>
<P>5,-4,1,1;
<p>
<p>
<P>-4,6,-4,2;
<p>
<p>
<P>1,-4,6,-3;
<p>
<p>
<P>
<p>
<p>
<P>输出:
<p>
<p>
<P>after the enhancedSquareRootMethod:the ans x rows is:(from x0 to xn-1)
<p>
<p>
<P>1.00000 1.00000 0.00000
<p>
<p>
<P>
<p>
<p>
<P>Test2:
<p>
<p>
<P>after the enhancedSquareRootMethod:the ans x rows is:(from x0 to xn-1)
<p>
<p>
<P>2.00000 1.00000 -1.00000
<p>
<p>
<P>
<p>
<p>
<P>网上关于这个主题的相关参考:
<p>
<p>
<P>http://jpkc.ecnu.edu.cn/gdds/xsxz/ZhangGengYun.htm
<p>
<p>
<P>http://www.ascc.net/pd-man/linpack/node14.html</P>
发表于 2006-4-15 19:57 | 显示全部楼层

回复:(tongji)是呀,文章中的变量说明也不全面

<DIV class=quote><B>以下是引用<I>tongji</I>在2006-4-14 16:04:39的发言:</B><BR>是呀,文章中的变量说明也不全面</DIV>
<br>具体的你可依找本数值算法的书看看<BR>比如丁丽娟的《数值计算方法》
发表于 2006-6-2 08:34 | 显示全部楼层

回复:(liljx_2008)[求助]求教大家,平方根辨识问题?...

<P><FONT color=#ff0000>liljx_2008、huhust各加威望1点,yejet加威望2点</FONT></P>
<P>多情清秋<BR>06.6.2</P>
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-27 06:57 , Processed in 0.127850 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表