线性方程组数值解法（2）

LU分解法：

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define n 6//n为矩阵维数

//矩阵三角分解A=LU

void Doolitte(double a[n][n],double L[n][n],double U[n][n])

{

int i,j,k;//分解L,U时,i,j,k在循环中计算循环次数，sum在循环中求和

double sum;

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

{

L[i][i]=1;

}

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

{

U[0][i]=a[0][i];//U的第一行元素

L[i][0]=a[i][0]/U[0][0];//L的第一列元素

}

//按公式计算L和U

for(k=1;k<n;k++)

{

for(i=k;i<n;i++)//计算U的第k行

{

sum=0;

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

{

sum=sum+L[k][j]*U[j][i];

}

U[k][i]=a[k][i]-sum;

if (U[k][k]==0)

{

printf("系数矩阵不满秩,无法进行LU分解!\n");

}

for(i=k+1;i<n;i++)//计算L的第k列

{

sum=0;

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

{

sum=sum+L[i][j]*U[j][k];

}

L[i][k]=(a[i][k]-sum)/U[k][k];

}

//输出L矩阵

printf("分解系数矩阵a得到的L矩阵为：\n");

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

{

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

{

printf("%-10f ",L[i][j]);

}

printf("\n");

}

printf("\n");

//输出U矩阵

printf("分解系数矩阵a得到的U矩阵为：\n");

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

{

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

{

printf("%-10f ",U[i][j]);

}

printf("\n");

}

printf("\n");

}

//L是一个下三角矩阵,且对角线元素为1,求解三角形方程组LY=b

void solve_L(double L[n][n],double y[n][1],double b[n][1])

{

int i,j,k;//j,k在循环中计算循环次数，n为矩阵维数，sum在循环中求和

double sum;

y[0][0]=b[0][0];

for(k=1;k<n;k++)

{

sum=0;

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

{

sum=sum+L[k][j]*y[j][0];

}

y[k][0]=b[k][0]-sum;

}

//U是一个上三角矩阵,求解三角形方程组UX=Y

void solve_U(double U[n][n],double x[n][1],double y[n][1])

{

int i,j,k;//j,k在循环中计算循环次数，n为矩阵维数，sum在循环中求和

double sum;

x[n-1][0]=y[n-1][0]/U[n-1][n-1];

for(k=n-2;k>=0;k--)

{

sum=0;

for(j=n-1;j>k;j--)

{

sum=sum+U[k][j]*x[j][0];

}

x[k][0]=(y[k][0]-sum)/U[k][k];

}

//输出解向量x

printf("解该方程组得到的x矩阵为：\n");

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

{

printf("%f\n",x[i][0]);

}

printf("\n");

}

//列向量归一化

void NormalMat(double x[n][1],double b[n][1])

{

int i;

double max;

max=fabs(x[0][0]);

for (i=1;i<n;i++)

{

if (max<fabs(x[i][0]))

{

max=fabs(x[i][0]);

}

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

{

b[i][0]=x[i][0]/max;

}

//输出迭代得到的向量b

printf("迭代得到的b向量为：\n");

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

{

printf("%-f\n",b[i][0]);

}

printf("\n");

}

int main()

{

double a[n][n]={{1,2,4,7,11,16},{2,3,5,8,12,17},{4,5,6,9,13,18},{7,8,9,10,14,19},{11,12,13,14,15,20},{16,17,18,19,20,21}};//输入系数矩阵

double b[n][1]={4,3,4,2,1,1};//自定义b1

double L[n][n]={0};//定义L矩阵

double U[n][n]={0};//定义U矩阵

double y[n][1];//定义y向量

double x[n][1];//定义x向量

int i,j,k;//i,j,k在循环中计算循环次数，n为矩阵维数，sum在循环中求和

double sum;

Doolitte(a,L,U);

solve_L(L,y,b);

solve_U(U,x,y);

for (i=1;i<10;i++)

{

printf("\n迭代%d次:\n",i);

NormalMat(x,b);//向量归一化

solve_L(L,y,b);//解方程组Ly=b

solve_U(U,x,y);//解方程组Ux=y

}

printf("\n\n结论:无论以任何初值出发，最终都将收敛到系数矩阵最大特征值(绝对值最大)对应的特征向量!\n");

return 0;

}

线性方程组数值解法（2）相关推荐

matlab算线性方程解,MATLAB计算方法3解线性方程组计算解法.pptx
第三章线性方程组数值解法解线性方程组 §3.1 直接法一. Gauss 消去法设有消元: 用Matlab实现顺序Gauss消去法在Matlab程序编辑器中输入:function x=nagauss ...
数值计算方法第三章—线性方程组的数值解法知识点总结
线性方程组的数值解法本文参考书为马东升著<数值计算方法> 高斯消去法顺序高斯消去法通过初等变换消去方程组系数矩阵主对角线以下的元素,而使方程组化为等价的上三角形方程组列主元高斯消去 ...
【计算方法】线性方程组的数值解法
题外话: 我上学期做了笔记的科目好像都炸了..像Java还有数据结构.计算方法也来做一个吧,反正迟早是要炸的一.综述线性方程组的解法可以分为两类:直接法和迭代法. 直接法通过有限四次运算得到精确解 ...
线性代数方程组数值解法
线性代数方程组的数值解法线性代数方程组数值解法一.向量范数与矩阵范数 1.1 向量范数 1.1.1 满足三个条件(向量范数公理) 1.1.2 常用的向量范数 1.2 矩阵范数 1.2.1 矩阵范数 ...
计算方法（六）：常微分方程初值问题的数值解法
文章目录常微分方程初值问题的数值解法欧拉(Euler)方法与改进欧拉方法欧拉方法欧拉公式的局部截断误差与精度分析改进欧拉方法龙格-库塔(Runge-Kutta)法构造原理经典龙格-库塔 ...
用matlab求解线性代数方程组,线性代数方程组数值解法与MATLAB实现综述
线性代数方程组数值解法及MATLAB 实现综述廖淑芳 20122090 数计学院 12计算机科学与技术1班(职教本科) 一.分析课题随着科学技术的发展,提出了大量复杂的数值计算问题,在建立电子计算 ...
通用求根算法zeroin_Modern Robotics运动学数值解法及SVD算法（C matlab）
前言原著之前CSDN已经注销,新CSDN Galaxy_Robot的博客_CSDN博客-机器人,C语言,我是谁?领域博主blog.csdn.net 这半个月的业余时间研究了机器人逆运动学的解析解法 ...
【数理知识】《数值分析》李庆扬老师-第9章-常微分方程初值问题数值解法
第8章回到目录无第9章-常微分方程初值问题数值解法 9.1 引言利普希茨 (Lipschitz) 条件 / 利普希茨常数定理1 解的存在唯一性定理定理2 解对初值依赖的敏感性 9.2 简单 ...
【数理知识】《数值分析》李庆扬老师-第7章-非线性方程与方程组的数值解法
第6章回到目录第8章第7章-非线性方程与方程组的数值解法 7.1 方程求根与二分法 7.2 不动点迭代法及其收敛性 7.3 迭代收敛的加速方法 7.4 牛顿法 7.5 弦截法与抛物线法 7.6 ...

线性方程组数值解法（2）

线性方程组数值解法（2）相关推荐

最新文章

热门文章