矩阵解线性方程组的程序

  cheney

在求解sin函数的样条差值的时候遇到需要解线性方程组,这个当然想到了用矩阵来求解。解矩阵的方法有很多。用程序解矩阵的想法一般有克莱姆法则、无回代的主元消去法、有回代的主元消去法以及其他方法。

1.三阶以下的矩阵一般用克莱姆法则来解

克莱姆法则表述如下,对于线性方程组:

wps_clip_image-23921

如果该方程组的系数行列式不等于零,即:

wps_clip_image-16079

那么该方程有唯一一组解:

wps_clip_image-9751wps_clip_image-5319wps_clip_image-18185wps_clip_image-12672wps_clip_image-14805。其中 wps_clip_image-10067是把系数行列式中的第 i 列元素用方程右端的常数项替换后得到的n阶行列式。

据此写出处理函数如下


	//求三阶行列式的
	double  Determinant3 (  double  list[ 3 ][ 3 ] )
	{
		double  sum1;
		double  sum2;
		
		sum1 = list[ 0 ][ 0 ]* list[ 1 ][ 1 ]*list[ 2 ][ 2 ];
		sum1 += list[ 0 ][ 1 ]* list[ 1 ][ 2 ]*list[ 2 ][ 0 ];
		sum1 += list[ 0 ][ 2 ]* list[ 1 ][ 0 ]*list[ 2 ][ 1 ];
		
		sum2 = list[ 0 ][ 2 ]* list[ 1 ][ 1 ]*list[ 2 ][ 0 ];
		sum2 += list[ 1 ][ 2 ]* list[ 2 ][ 1 ]*list[ 0 ][ 0 ];
		sum2 += list[ 2 ][ 2 ]* list[ 0 ][ 1 ]*list[ 1 ][ 0 ];

		return  (sum1-sum2);
	}
	
	void  Matrix3(  double  buff_A[ 3 ][ 3 ] ,  double  buff_B[ 3 ],  double  s[ 3 ] )
	{
		double  buff_C[ 3 ][ 3 ];
		double  v;
		int  j;
		memcpy( buff_C,buff_A, sizeof ( double ) *9 );
		
		// 系数行列式的值
		v = Determinant3(buff_C);

		if ( v <  0 . 0000000000000001  );
		printf("系数不能为零");

		// 替换第一列
		for (j= 0  ; j < 3 ;j++ )
		{
			buff_C[j][ 0 ]=buff_B[j];       
		}
		s[ 0 ]= Determinant3(buff_C)/v;

		memcpy(buff_C,buff_A, sizeof ( double ) *9 );
		
		// 替换第二列
		for (j= 0  ; j < 3 ;j++ ){
		buff_C[j][ 1 ]=buff_B[j];       
		}
		s[ 1 ]= Determinant3(buff_C)/v;

		memcpy(buff_C,buff_A, sizeof ( double ) *9 );
		
		// 替换第三列
		for (j= 0  ; j < 3 ;j++){
		buff_C[j][ 2 ]=buff_B[j];       
		}
		s[ 2 ]= Determinant3(buff_C)/v;

2.消元法

对于三阶以上的矩阵,用行列式法就不太简单了,因为四阶以上行列式的要成的项越来越多,此时,一般还是选用还是消元法。

对于线性方程组:

wps_clip_image-29267

可以列出增广矩阵:

wps_clip_image-8727

然后用消元法使系数矩阵转化为单位矩阵:

wps_clip_image-16109

则右边的系数变为方程的解。

基于这种思想的用程序如下


	// 寻找主元
	void  find ( int  n,  double  *buff_X , int  i )
	{
		double  c = fabs( *(buff_X + i*(n + 1 )+i ) );
		int  k = i;
		int  j;
		for (j= i+ 1 ;j < n ; j++)
		{
			if (   fabs( *(buff_X + j*(n+ 1 ) + i)  > c )
			{
				c = fabs( *(buff_X + j*(n+ 1 ) + i );
			}
		}  

		k = j;

		if (k != i)	// 交换行
		{
			for ( j= 0  ; j <= n ; j++ )
			{
				c= *(buff_X + k*(n+ 1 ) + j) ;
				*(buff_X + k *(n+ 1 ) + j) = *(buff_X+ i*(n+ 1 ) + j);
				*(buff_X + i*(n+ 1 ) + j) = c;
			}
		}
	}
	
	// 各行的系数除以主元
	void  div( int  n, double  *buff_X,  int  i)
	{
		int  j;
		double  c;
		c = *(buff_X + i*(n+ 1 ) + i );
		*(buff_X + i*(n+ 1 ) + i )= 1 . 0 
		
		for (j= i + 1  ; j <= n ; j++)
		{
			*(buff_X + i*(n+ 1 ) +  j) /=c ;
		}
	}
	
	// 消元
	void  del( int  n , double  *buff_X , int  i )
	{
	int  j,k;
	float  c;
	for (j=i+ 1  ; j < n; j++ );
	{               
		if (*(buff_X + j*(n+ 1 ) + i )
		c = *(buff_X + j*(n+ 1 ) + i )
		*(buff_X + j*(n+ 1 )  +i ) = 0 
		for ( k=i+ 1  ;k&lt;=n ;k++ )
		{
			*(buff_X  +j*(n+ 1 )  +k ) -= c* (*(buff_X +i*(n+ 1 ) +k )
		}
	}	

void  Matrix( int  n,  double  *buff_X,  double  *s )
{
	int  j;
	int  i;
	for (j = 0  ; j < n ; j++ )
	{
		//对于每一列
		//寻找主元
		if (j < n - 1 )
		{
			find(n,buff_X,j); 
		}

		div(n,buff_X,j)
		
		if (j < n - 1 )
		{ 
			del(n,buff_X,j);    
		}
	}
	// 回代
	for (i = n- 2 ;i <= 0 ; i-- )
	{
		for (j = n- 1 ;j < i; j--)
		{
			*( buff_X + i*(n+ 1 ) + n ) -= *( buff_X + i*(n+ 1 ) + j) * (*(buff_X + j*(n+ 1 ) + n));
		}
	}
	
	for ( i = 0  ; i < n ; i++ )
	{ 
		*(s+i) = *(buff_X + i*(n+ 1 ) + n )
	}
   
}

3.调用Matlab解矩阵

还有一种比较懒的方法解矩阵,那就是调用Matlab后台计算,这个方法不能用在嵌入式系统里了。以下是VS2008连接Matlab 7.0的操作。都不是新版本,即使不同版本也大同小异。

3.1配置

通过菜单工具/选项,打开选项页,点击项目和解决方案,然后在页面右面“显示以下内容的目录”下拉列表框中选择“包含文件”,添加路径:“d:/Program files/MATLABR2007a/extern/include/”。

选择“库文件”,添加路径:“d:/Program Files/MATLABR2007a/extern/lib/win32/microsoft/”。具体路径可能不同,主要要找到某个版本下的 extern 文件夹。

右击工程/属性,打开项目属性页,选择链接器/输入,在附加依赖项编辑框中,添加文件名 libmx.lib、libmat.lib、libeng.lib。

在工程中添加#include "..matrix_c.h"要给出具体路径。

####3.2 常用 API

//在调用Matlab引擎之前,首先应在相关文件中加入一行:#include "enging.h",该文件包含了引擎API函数的说明和所需数据结构的定义

//Matlab7.0 中可以在VC中调用的引擎函数分别如下

	/*
	启动Matlab进
	参数startcmd是用来启动Matlab引擎的字符串参数,在Windows操作系统中只能为NULL
	函数返回值是一个Engine类型的指针,它是在engine.h中定义的engine数据结构
	*/
	
	Engine *engOpen (
		const   char  *startcm
	)
	
	
	/*
	关闭Matlab 引
	参数ep代表要被关闭的引擎指针
	*/
	
	int  engClose(
		Engine    *ep         /* engine pointer *
	)
	
	/*
	发送命令让Matlab执
	参数ep为函数engOpen返回的引擎指针,字符串string为要matlab执行的命令
	函数返回值为0表示成功执行,返回1说明执行失败(如命令不能被Matlab正确解释或Matlab引擎已经关闭了)
	*/
	
	int  engEvalString
		Engine    *ep,       
		const   char  *strin
	)
	
	
	
	/*
	参数ep为Matlab引擎指针,p为用来保存输出结构的缓冲区,n为最大保存的字符个数,通常就是缓冲区p的大小
	该函数执行后,接下来的engString函数所引起的命令行输出结果会在缓冲区p中保存
	*/
	
	int  engOutputBuffer
		Engine    *ep,   
		char     *buffer,   
		int      bufle
	)
	
	
	/*
	从Matlab引擎工作空间中获取变量
	参数ep为打开的Matlab引擎指针,name为以字符串形式指定的数组名
	函数返回值是指向name数组的指针,类型为mxArray*
	*/
	
	mxArray *engGetVariable
		Engine    *ep,   
		const   char  *name   
	)
	
	
	/*
	向Matlab引擎工作空间写入变量
	参数ep为打开的Matlab引擎指针,mp为指向被写入变量的指针,name为变量写入后在Matlab引擎工作空间中的变量名
	函数返回值为0表示写入变量成功,返回值为1表示发生错误
	*/
	
	int  engPutVariable
		Engine    *ep,      
		const   char  *var_name
		const  mxArray *ap  
	)
	
	
	/*
	以engine方式调用Matlab的时候,会打开Matlab主窗口,可在其中随意操作。但有时也会干扰应用程序的运行,可用以下设置是否显示该窗口
	参数ep为打开的Matlab引擎指针,value为是否显示的标志,取值true(或1)表示显示Matlab窗口,取值false(或0)表示隐藏Matlab窗口
	函数返回值为0表示设置成功,为1表示有错误发生
	*/
	
	int  engSetVisible(
		Engine *ep,      
		bool  newVal
	)
	
	/*
	要获得当前Matlab窗口的显示/隐藏情况,可以调用函数
	参数ep为打开的Matlab引擎指针,Value为用来保存显示/隐藏情况的变量(采用指针方式传递)
	函数返回值为0表示获取成功,为1表示有错误发生
	*/
	
	int  engGetVisible(
		Engine *ep,     
		bool * bVal
	)
	
	
	/*
	创建一个单次使用的Matlab进
	参数startcmd是用来启动Matlab引擎的字符串参数,在Windows操作系统中只能为NULL
	参数reserved是为方便以后扩展而设置,只能是NULL
	参数retstatus为状
	*/
	
	Engine *engOpenSingleUse(
		const   char  *startcmd,
		void  *reserved
		int  *retstatus
	)

3.3示例程序

为了以后使用方便,将API再次封装,matrix_c.cpp

// 增广矩阵求解
// ep 为有效的引擎指
// n 为系数矩阵阶
// buff_A 为系数矩阵指
// buff_B 为常数矩阵指
// S 指向解的空
return 执行正确返回
int  SolMatrix(Engine *ep,  int  n,  double  *buff_A,  double  *buff_B,  double  *S)
{
	int  err= 0 ;
	//定义一个n行n列的实数系数
	mxArray *A = mxCreateDoubleMatrix(n,n, mxREAL);
	//定义一个1行n列的实数常数
	mxArray *B = mxCreateDoubleMatrix( 1 ,n, mxREAL);
	//结果的缓
	mxArray *Sol = mxCreateDoubleMatrix( 1 ,n, mxREAL);
	//拷贝数据进入结
	memcpy(mxGetPr(A), buff_A , n*n *sizeof ( double ));
	memcpy(mxGetPr(B), buff_B , n *sizeof ( double ));
	
	//将数据输入MatLa
	err += engPutVariable(ep, "A",A);
	err += engPutVariable(ep, "B",B);
	//计算
	err += engEvalString(ep, "S=B/A");
	//获取结果
	Sol = engGetVariable( ep,"S");
	if  ( NULL  == Sol )
	{
		err++;
	}   
	if ( 0 ==err )
	{
		memcpy( S, mxGetPr(Sol) , n *sizeof ( double ));
		return   0;
	}
	else
	{
		return  err;
	}
}




marix_c.h

/**
文件:矩阵求
作者:Channe
时间:2012-06-2
版本:V 0.
说明:如果安装了Matlab,开启定义__MATLAB(默认)。如果没有安装Matlab,则使自定义的解矩阵的方法
*/

#ifndef __MATRIX_C_H_

#define __MATRIX_C_H_


//默认设

//如果未安装Matlab 注释掉此

#define __MATLA


#ifdef __MATLA

#include "engine.h

#else

#define __EMBED

#endif


int  SolMatrix( int  n,  double  *buff_X,  double  *s);

int  SolMatrix(Engine    *ep,  int  n,  double  *buff_A,  double  *buff_B,  double  *S);

#endif

main.cpp

/**
工程:实现VS2008 和 Matlab 的连
作者:Channe
时间:2012-06-2
*/

#include <iostream>
#include <stdio.h>
#include <math.h>
#include "matrix_c.h"
using   namespace  std;


void  main()
{
	const   int  N =  50 ;
	double  x[N],y[N];
	int  j =  1 ;
	for  ( int  i= 0 ; i < N; i++) //计算数组x和y
	
	x[i] = ( i + 1 )
	y[i] = sin(x[i]) + j * log(x[i])
	j *= - 1 

	//定义Matlab引擎指针
	Engine *ep;
	if  (!(ep=engOpen( NULL ))) //测试是否启动Matlab引擎成功
	{
		exit( 1 );
	}

	double  buff_A[ 2 ][ 2 ] = { 2 , 1 , 1 , 2 };
	double  buff_B[ 2 ]= { 1 , 3 };
	double  s[ 2 ]= { 0 , 0 };
	SolMatrix( ep ,  2  , ( double  *)buff_A, ( double  *)buff_B, ( double  *)s);
	cout << s[ 0 ] << "," << s[ 1 ] << endl;
	printf("%.15f,%.15fn",s[ 0 ],s[ 1 ]);
	cout << "执行完毕" <<endl;
	cin.get();
	//关闭Matlab引擎
	engClose(ep);
}

## 扩展阅读

《 精通Matlab与C混合编程》