精华内容
下载资源
问答
  • 矩阵列主元三角分解

    2020-12-18 21:57:38
    一,列主元三角分解定理 如果A为非奇异矩阵,则存在排列矩阵P使 PA=LU,其中L为下三角矩阵,U为上三角矩阵,即A=P-1LU 二,列主元三角分解Python代码 # 自己原创 def pivot_lu_decomposition(coefficient_matrix: np....

    一,列主元三角分解定理

    如果A为非奇异矩阵,则存在排列矩阵P使 PA=LU,其中L为下三角矩阵,U为上三角矩阵,即A=P-1LU

    二,列主元三角分解Python代码

    # 自己原创
    def pivot_lu_decomposition(coefficient_matrix: np.ndarray, right_hand_side_vector: np.ndarray):
        """
        实现方程Ax=b系数矩阵A的pivoted LU decomposition
        :param coefficient_matrix: 初始系数矩阵A
        :param right_hand_side_vector: 初始常数列向量b
        :return: 排列矩阵P,单位下三角矩阵L,上三角矩阵U,常数列向量b
        """
        # first step: evaluate the lu decomposition condition
        rows, columns = coefficient_matrix.shape
        if rows == columns:  # judge if it is a square matrix
            for k in range(rows):  # 判断各阶顺序主子式是否为0
                if det(coefficient_matrix[:k + 1, :k + 1]) == 0:
                    raise Exception("cannot generate LU decomposition")
            else:
                # pivot LU decomposition
                lower_triangular_matrix = np.eye(rows)
                permutation_matrix = np.eye(rows)
                for k in range(rows - 1):
                    max_index = np.argmax(abs(coefficient_matrix[k:, k]))
                    coefficient_matrix[[k, max_index + k], :] = coefficient_matrix[[max_index + k, k], :]
                    right_hand_side_vector[[k, max_index + k], :] = right_hand_side_vector[[max_index + k, k], :]
                    permutation_matrix[[k, max_index + k], :] = permutation_matrix[[max_index + k, k], :]
                    for row in range(k + 1, rows):
                        multiplier = coefficient_matrix[row, k] / coefficient_matrix[k, k]
                        coefficient_matrix[row, k:] += -multiplier * (coefficient_matrix[k, k:])
                        right_hand_side_vector[row] += -multiplier * right_hand_side_vector[k]
                        lower_triangular_matrix[row, k] = multiplier
        else:
            raise Exception("ERROR:please pass a square matrix.")
    
        return permutation_matrix, lower_triangular_matrix, coefficient_matrix, right_hand_side_vector
    
    展开全文
  • 列主元三角分解法在matlab中的实现.doc
  • LU分解法,列主元三角分解法MATLAB代码,有详细注释,易于理解
  • LU分解法且是列主元三角分解法MATLAB代码,有详细注释,顺着思路看容易理解
  • 最近学了一般线性方程组的直接解法——列主元三角分解。于是我写了一个程序,用来获取三角分解后的上三角矩阵、下三角矩阵以及方程的解。 当然这些用matlab去实现更简单,但是用C的话更能够感受到其中的算法流程。 ...

    最近学了一般线性方程组的直接解法——列主元三角分解。于是我写了一个程序,用来获取三角分解后的上三角矩阵、下三角矩阵以及方程的解。

    当然这些用matlab去实现更简单,但是用C的话更能够感受到其中的算法流程。

    由于C语言进行矩阵元素的话比较麻烦,所以我把它写在一个叫array_number.txt文件里了。

     方便读取~

    代码:

    #include<stdio.h>
    #include<stdlib.h>
    #include<math.h>
    double  **matrix_input(int n,FILE *file_pointer)//从指针file_pointer指向的文件里边读入二维数组信息 
     {
     	int i,j;
    	double  **array;
    	array=(double **)calloc(n,sizeof(double *));
    	for(i=0;i<n;i++)
    	{
    	    array[i]=(double *)calloc(n+1,sizeof(double));//增广矩阵 
    	    for(j=0;j<n+1;j++)
    		fscanf(file_pointer,"%lf",&array[i][j]);
        }
    	return array; 
     }
     void matrix_printf(double**array,int n) //写一个矩阵显示函数,方便观察矩阵元素的变化 
    {
    	int i,j;
    	for(i=0;i<n;i++)
    	{
    		for(j=0;j<n+1;j++)
    		{
    		    printf(" %lf ",array[i][j]);
    		}
    		 printf("\n");
    	}
    }
    void Guass(double**array_U,int n) //改进的gauss消去  三角分解   
    {
    	int i,j,k; 
    	int *p;
    	double m;
        double **array_L;
    	p=(int *)calloc(n,sizeof(int));
    	array_L=(double **)calloc(n,sizeof(double *));//声明一个二维动态数组用于存储下三角矩阵 
    	for(i=0;i<n;i++)
    	{
    	    array_L[i]=(double *)calloc(n+1,sizeof(double));
    	    for(j=0;j<n;j++)array_L[i][j]=0;
        }
        for(i=0;i<n;i++)//下三角矩阵的对角线元素为1 
        {
        	array_L[i][i]=1;
    	}
    	for(i=0;i<n-1;i++)  //进行列主元三角分解 
    	{
    		p[i]=i; 
    		m=array_U[i][i];
    		//选主元
    		for(j=i+1;j<n;j++)
    		{
    			if(fabs(m)<fabs(array_U[j][i])) //寻找该列之中最大的 
    			{
    				m=array_U[j][i];
    				p[i]=j;
    				printf("该列最大元素:%lf\n",m);
    			}
    		} 
    	    
    		if(p[i]!=i)//考虑换行 
    		{
    			for(j=0;j<n+1;j++)
    			{
    				m=array_U[i][j];
    				array_U[i][j]=array_U[p[i]][j];
    				array_U[p[i]][j]=m;
    			}
    			printf("换行后:\n"); 
    		    matrix_printf(array_U,n); 
    		}
    		//消元
    		for(j=i+1;j<n;j++)
    		{
    			array_U[j][i]/=array_U[i][i];
    			array_L[j][i]=array_U[j][i];   //将存储在array_U中的下三角矩阵数据拿出 
    			for(k=i+1;k<n+1;k++)
    			array_U[j][k]-=array_U[j][i]*array_U[i][k];
    			array_U[j][i]=0;
    			printf("消元后:\n");
    			matrix_printf(array_U,n); 
    		}
    	}
    	if(array_U[n-1][n-1]==0)
    	{
    		printf("该矩阵有问题,可能为奇异矩阵\n");
    		exit(0);
    	}
    	printf("三角分解后下三角矩阵L:\n");
    	matrix_printf(array_L,n); 
    }
    
    double *changearrayb(double**array_U,int n)//将增广矩阵的最后一列存储到数组b中 
    {
    	int i,j;
    	double *b;
    	b=(double *)calloc(n,sizeof(double));
    	for(i=0;i<n;i++)
    	{
    	    b[i]=array_U[i][n];	
    	}
    	return b;
    }
    double UpTriangle(double** L, double* b, int n) //上三角矩阵求解  回代法 
    {
    	int i, j;
    	for (i = n - 1; i >= 0; i--)
    	{
    		for (j = i + 1; j <= n - 1; j++)
    		{
    			b[i] -= L[i][j] * b[j];
    		}
    		b[i] = b[i] / L[i][i];
    	}
    	return *b;
    }
    int main()
    {
    	int n,i,j;
    	double **array_U,*b;
    	FILE *fp=fopen("array_number.txt","r");
    	fscanf(fp,"%d",&n);
    	array_U=matrix_input(n,fp);//已读入
    	Guass(array_U,n);
    	b=changearrayb(array_U,n);
    	*b=UpTriangle(array_U,b,n);
    	for(i=0;i<n;i++)
    	{
    		printf("X%d:%lf\n",i+1,b[i]);
    	}
    	system("pause");
    } 
    

    这其中用了大量的二级指针和二维数组,存在以指针方式传递矩阵数据的方法。我写完之后感觉对指针、数组以及以指针和数组为参数的函数调用理解更深刻了。

    我使用的是VC6.0,如果使用其他的编译器,我不知道会不会有差异。

    展开全文
  • 用列主元高斯消去法和列主元三角分解法解线性方程 【目的和意义】 高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。 用...

    用列主元高斯消去法和列主元三角分解法解线性方程

    【目的和意义】

    高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。

    用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵(上三角矩阵、单位矩阵等),而三角形方程组则可以直接回带求解

    用高斯消去法解线性方程组
    除法运算步骤为

    ,加减运算步骤为

    。相比之下,传统的克莱姆法则则较为繁琐,如求解20阶线性方程组,克莱姆法则大约要

    次乘法,而用高斯消去法只需要3060次乘除法。

    在高斯消去法运算的过程中,如果出现abs(A(i,i))等于零或过小的情况,则会导致矩阵元素数量级严重增长和舍入误差的扩散,使得最后的计算结果不可靠,所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法,从而使计算结果更加精确。

    2、列主元三角分解法

    高斯消去法的消去过程,实质上是将A分解为两个三角矩阵的乘积A=LU,并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中,矩阵L和U可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法

    采用选主元的方式与列主元高斯消去法一样,也是为了避免除数过小,从而保证了计算的精

    disp(‘原方程为AX=b:’) %显示方程组

    A
    b
    disp(’------------------------’)
    n=length(A);
    for k

    展开全文
  • [数值算法]列主元三角分解

    千次阅读 2005-09-14 12:12:00
    [数值算法]列主元三角分解法 By EmilMatthew 05/9/13 列主元三角分角法是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,就是避免小数作为分母项. 有

    [数值算法]列主元三角分解法

                                                                                                      By EmilMatthew 05/9/13

           列主元三角分角法是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,就是避免小数作为分母项.

           有一些麻烦是的是要对原矩阵的增广矩阵部分进行LU分解,而且每次要对当前未处不景行作求列主行的下标,然后做变换.

    其它的注意点和LU分解法其本一致,而最后求出U矩阵后,可以直接用回代法求出解x,而不必先求y再用Ux=y来求解向量x.

    其实,这里根本没有必要用到数学上推导时用的LU两个三角矩阵,只要用一个增广矩阵就可以解决战斗了,不过我实现时还是用了老方法,以后有待改进.

    /*

    LUColMainMethod, coded by EmilMathew 05/9/13, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.

    */

     

    void LUColMainMethod(Type** matrixArr,Type* bList,Type* xAnsList,int len)

    {

           /*Core think of this algorithm:

           matrix_L*yAnsList=bList;

           matrix_U*xAnsList=yAnsList;

           */

           Type** matrix_L,** matrix_U;/*The core two matrix object of the LU method.*/

           Type* yAnsList;

          

           int i,j,k;/*iterator num*/

           int i2,j2;

           Type sum;/*the temp tween data*/

           /*LU Col Main Spec Data*/

           int maxRowIndex;

           Type** matrixA_B;

     

     

           /*input assertion*/

           assertF(matrixArr!=NULL,"in LUPationMethod,matrixArr is NULL/n");

           assertF(bList!=NULL,"in LUPationMethod,bList is NULL/n");

           assertF(xAnsList!=NULL,"in LUPationMethod,xAnsList is NULL/n");      

          

           /*memory apply*/

           matrix_L=(Type**)malloc(sizeof(Type*)*len);

           twoDArrMemApply(&matrix_U,len,len+1);

          

           for(i=0;i<len;i++)

                         {

                                matrix_L[i]=(Type*)malloc(sizeof(Type)*len);

                         //     matrix_U[i]=(Type*)malloc(sizeof(Type)*len);

                         }

          

           yAnsList=(Type*)malloc(sizeof(Type)*len);

           /*==end of memory apply==*/

          

           /*---assertion after apply the mem---*/

           assertF(matrix_L!=NULL,"in LUPationMethod,matrix_L is NULL/n");

           assertF(matrix_U!=NULL,"in LUPationMethod,matrix_U is NULL/n");            

           assertF(yAnsList!=NULL,"in LUPationMethod,yAnsList is NULL/n");

           /*---end of assertion---*/

    //     printf("arr mem applyed/n"); 

           /*----data initialization----*/

           for(i=0;i<len;i++)

           {

                  for(j=0;j<len;j++)

                  {

                         matrix_L[i][j]=0;

                         matrix_U[i][j]=0;

                  }

                  matrix_U[i][j]=0;

           }

     

           for(i=0;i<len;i++)

                  matrix_L[i][i]=1;

           /*matrix A_B prepare*/

           twoDArrMemApply(&matrixA_B,len,len+1);

     

           matrixCopy(matrixArr,matrixA_B,len,len);

           for(i=0;i<len;i++)

                  matrixA_B[i][len]=bList[i];

           printf("show copied matrix:/n");   

           show2DArrFloat(matrixA_B,len,len+1);

          

    for(i=0;i<len;i++)

                  {

                         maxRowIndex=getMaxRowIndexLU(matrixA_B,i,matrix_L,matrix_U,len);

                 

                         if(i!=maxRowIndex)

                         {

                                swapTwoRow(matrixA_B,i,maxRowIndex,len+1);

                                swapTwoRow(matrix_L,i,maxRowIndex,len);

                                swapTwoRow(matrix_U,i,maxRowIndex,len);

                         }

     

                         /*matrixU make*/

                         for(j=i;j<len;j++)

                         {

                                sum=0;

                                for(k=0;k<i;k++)

                                       sum+=matrix_L[i][k]*matrix_U[k][j];

                                matrix_U[i][j]=matrixA_B[i][j]-sum;

                         }

                        

                         matrix_U[i][j]=matrixA_B[i][j]-sumArr1_IKByArr2_KJ(matrix_L,matrix_U,i,j,0,i-1);

                         matrixA_B[i][j]=matrix_U[i][j];

                        

                         /*matrixL make*/

                         if(i<len-1)

                         {

                                j2=i;

                                for(i2=j2+1;i2<len;i2++)

                                {

                                       sum=0;

                                       for(k=0;k<j2;k++)

                                              sum+=matrix_L[i2][k]*matrix_U[k][j2];

                                       matrix_L[i2][j2]=(matrixA_B[i2][j2]-sum)/matrix_U[j2][j2];

                                }

                         }

                        

                        

                  }

           //copyBack

     

           for(i=0;i<len;i++)

                  bList[i]=matrixA_B[i][len];

           printf("show bList/n");

           showArrListFloat(bList,0,len);

           /*Adjusting of matrix_L*/

           for(i=0;i<len;i++)

           {

                  for(j=i;j<len;j++)

                  {

                         if(i==j)

                                matrix_L[i][i]=1;

                         else

                                matrix_L[i][j]=0;

                  }

           }

     

           printf("matrix L/n");

           show2DArrFloat(matrix_L,len,len);

           printf("matrix U/n");

           show2DArrFloat(matrix_U,len,len);

     

           for(i=len-1;i>=0;i--)

           {

                  xAnsList[i]=(bList[i]-sumArr_JKByList_K(matrix_U,xAnsList,i,i+1,len-1))/matrix_U[i][i];

           }

           /*                 

           downTriangleMethod(matrix_L,bList,yAnsList,len);

          

           upTriangleMethod(matrix_U,yAnsList,xAnsList,len);

           */

     

           /*memory free*/

          

           for(i=0;i<len;i++)

                  {

                         free(matrix_L[i]);

                         free(matrix_U[i]);

                  }

                 

           free(matrix_L);

           free(matrix_U);

           free(yAnsList);

           free(matrixA_B);

    }

          

    //辅助函数,求列主元的行位置.

    int getMaxRowIndexLU(Type** inMatrixArr,int currentI,Type** matrix_L,Type** matrix_U,int size)

    {

           int i,maxRowIndex;

           Type tmpSMax,tmpCompare;

                               

           assertF(inMatrixArr!=NULL,"in getMaxRowIndexLU,inMatrixArr is NULL/n");

           assertF(matrix_L!=NULL,"in getMaxRowIndexLU,matrix_L is NULL/n");

           assertF(matrix_U!=NULL,"in getMaxRowIndexLU,matrix_U is NULL/n");

                               

           tmpSMax=0;

    //     maxRowIndex=currentI;

           for(i=currentI;i<size;i++)

           {

                  tmpCompare=(float)fabs(inMatrixArr[i][currentI]-sumArr1_IKByArr2_KJ(matrix_L,matrix_U,i,currentI,0,currentI-1));

                  if(tmpCompare>tmpSMax)

                  {

                         tmpSMax=tmpCompare;

                         maxRowIndex=i;

                  }

           }

          

           return maxRowIndex;

    }

     

    我所写的程序的源码下载:

    http://emilmatthew.51.net/downloads/luColMain.rar

    //如果上面这个链接无法响应下载(有可能是被网站给屏蔽掉了),则可使用下载工具(如迅雷等)下载

     

    测试结果:

    test1:

    input:

    3;

    4,-2,-4,10;

    -2,17,10,3;

    -4,10,9,-7;

     

    output:

    after the LUPationMethod:the ans x rows is:(from x0 to xn-1)

       2.00000    1.00000   -1.00000

     

    test2:

    input:

    3;

    12,-3,3,15;

    18,-3,1,15;

    -1,2,1,6;

     

    after the LUPationMethod:the ans x rows is:(from x0 to xn-1)

       1.00000    2.00000    3.00000

    展开全文
  • 直接利用矩阵的列主元三角分解求解线性方程组Ax=b的计算过程可按如下步骤进行:第一步:运用算法1.2.2计算A的列主元LU分解:PA=LU;第二步:运用算法1.1.1解下三角形方程组Ly=Pb;第三步:运用算法1.1.2解上三角形...
  • 推荐答案上善若水666来自团队: ...这个定理说明先对A进行对换矩阵的行得到PA,然后再对PA进行LU分解是可行的.证明如下:A选主元的LU分解实际是对应这样的矩阵相乘U=(Ln-1En-1)..(L2E2)(L1E1)A看等号右边我们来解释一...
  • (1)编制解n阶线性方程组Ax=b的列主元Gauss消去法的通用程序;
  • 例如,我们知道上三角矩阵和下三角矩阵是容易求解的,或者对角矩阵是最理想的求解形式,我们通过矩阵的分解去将原本的复杂问题,转化为若干个易于求解的矩阵来运算。一、高斯消去法与LU分解大学本科期间学校讲授的...
  • Gauss消去法 列主元Gauss消去法 ...特别的,已知线性方程组,其中矩阵A为三对角阵,应用三角分解,可得 周期三对角方程组: 转载于:https://www.cnblogs.com/wander-clouds/p/9991481.html...
  • Doolittle三角分解法,与上面的Guass列选主元消去法有着一些联系,至少前面的输入代码是差不多的:题目: 2、Doolittle三角分解法 [ 2 10 0 -3 ] [-3 -4 -12 13] [ 1 2 3 -4 ] [ 4 14 9 -13]
  • Doolittle分解法(三角分解算法)求解线性方程组(MATLAB实现) 求解线性方程组AX=b, b=(78,75,101 ,35 ,72,91 ,73 ,39 ,76,129) ,其中系数矩阵A如下: 1、三角分解算法 设线性方程组AX=B的系数A矩阵存在三角...
  • 本程序是用C++语言编写。本程序中LU分解采用的是列主元三角分解法,即PLU分解,实现LU=PA。也是中科院矩阵分析程序作业。
  • 这次是这周学的三角分解,matlab有对应的函数,在这里主要记录一下大概的分解结果和方法,最后自己用代码实现了一下。具体的求解方法和直接三角分解的方法一致,都是写出LU的形式,然后展开LU并与A对应,求取L和U的...
  • 列主元Gauss消去法(C++实现)

    千次阅读 2018-11-18 21:38:55
    目的:编写解n阶线性方程组AX=b的列主元三角分解法的通用程序; 原理:列主元素消去法是为控制舍入误差而提出来的一种算法,列主元素消去法计算基本上能控制舍入误差的影响,其基本思想是:在进行第 k(k=1,2,...,n-...
  • 《矩阵论》学习笔记(四)-1:4.1 矩阵的三角分解 矩阵的三角分解1.一般方阵- 高斯消去法与矩阵的LU分解2.可逆矩阵- Doolittle/Crout分解3.分块方阵- 拟LU分解与拟LDU分解 文章目录《矩阵论》学习笔记(四)-1:4.1 ...
  • 1. 矩阵三角分解法 实现代码 2. LU分解 2.1 基本步骤 2.2 LU分解的计算公式 2.3 LU分解的结果表示 实现代码 3. 选主元的LU分解 3.1 基本原理 3.2 算法过程描述 实现代码 4. 对称正定矩阵的Cholesky分解 4.1 基本原理...
  • Python实现列主元高斯消去法与LU分解法 数值分析:Python实现列主元高斯消去法与LU分解法求解线性方程组 ...#目的:熟悉列主元消去法,以及三角分解法等直接求解线性方程组的算法 #列主元消元法 def
  • 最近数值计算学了Guass列主消元法和三角分解法解线性方程组,具体原理如下: 1、Guass列选主元消去法对于AX=B 1)、消元过程:将(A|B)进行变换为,其中是上三角矩阵。即: k从1到n-1 a、列选主元 选取第k列...
  • function [L,U,y,x] = Doolittle(A, b) [row_a, col_a] = size(A); for j = 1:col_a U(1,j) = A(1,j); end L(1,1) = 1; for i = 2:row_a L(i,1) = A(i,1)/A(1,1); end for i = 2:row_a ... for k = 1:i-1
  • 计算方法实验(高斯消除法,三角分解法,最小二乘法,等等) 还有其他的,包挂报告,
  • 数值分析3-解线性方程组迭代法解线性代数方程组直接法高斯消元,高斯主元素消元矩阵三角分解法,列主元三角分解平方根法追赶法方程组性态和误差分析 解线性代数方程组直接法 高斯消元,高斯主元素消元 矩阵三角分解...
  • 一、实验内容二、代码(python)... 列主元高斯消元法 A:系数增广矩阵 n:未知数个数 ''' def main_element_gauss(A,n): for i in range(0,n-1): if(np.max(A[i:,i])!=A[i,i]): #如果当前系数不是最大值,则列主元
  • 每个代码都可以运行哦 运行环境我的是VC6.0
  • 列主元消去法和LU矩阵分解以及矩阵迭代求解

    多人点赞 热门讨论 2021-06-06 13:56:22
    三角求解 wei #上三角求解: def up_solve(A,b): n = len(A[0]) A[n-1,n-1] = b[n-1][0]/A[n-1,n-1] for i in range(1,n): k = n-i-1#取逆序 A[k,k] = (b[k][0] - sum(A[k]*A[k+1]))/A[k,k] A[k]=np....

空空如也

空空如也

1 2 3 4 5 ... 14
收藏数 268
精华内容 107
关键字:

列主元三角分解