# include<stdio.h>
# include<math.h>
#include<conio.h>
# define n 3
int main ()
{ int i,j,k,row,v,t;
int M[n];
float m,r,max0,max1;
float a[n][n],b[n],s[n],l[n][n],u[n][n],x[n],y[n];
a[0][0]=0.5,a[0][1]=1.1,a[0][2]=3.1;
a[1][0]=5.0,a[1][1]=0.96,a[1][2]=6.5;
a[2][0]=2.0,a[2][1]=4.5,a[2][2]=0.36;
b[0]=6.0,b[1]=0.96,b[2]=0.02;
for(k=0;k<n;k++)
{ for(i=k;i<n;i++)
{ m=0;
if(k-1>=0)
{ for(t=0;t<k;t++)
m=m+l[i][t]*u[t][k];
}
s[i]=a[i][k]-m;
}
max0=abs(s[k]);
for(i=k;i<n;i++)
{ if(abs(s[i])>max0)
{ max0=abs(s[i]);
row=i;}
}
### 选主元的Doolittle分解法程序详解
#### 一、引言
Doolittle分解法是一种求解线性方程组的有效方法，它将一个方阵分解为一个下三角矩阵\(L\)和一个上三角矩阵\(U\)的乘积。在实际应用中，为了提高计算的稳定性，通常会采用选主元策略来避免计算过程中出现的除零或误差累积问题。本文将详细解析一段实现选主元Doolittle分解法的C语言程序，并深入探讨其中的关键概念和技术细节。
#### 二、基础知识
1. **Doolittle分解**：对于一个\(n \times n\)的可逆矩阵\(A\)，存在唯一的下三角矩阵\(L\)（对角线元素均为1）和上三角矩阵\(U\)，使得\(A = LU\)。
2. **选主元技术**：在分解过程中选择绝对值最大的元素作为当前列的主元，以减少数值不稳定性的风险。
3. **求解过程**：
- **分解阶段**：首先进行矩阵分解得到\(L\)和\(U\)。
- **前向代换**：利用\(L\)矩阵和右侧向量\(b\)求得中间向量\(y\)。
- **后向代换**：再利用\(U\)矩阵和\(y\)向量求得最终解向量\(x\)。
#### 三、代码解析
##### 1. 初始化与定义
```c
# include<stdio.h>
# include<math.h>
#include<conio.h>
# define n 3
```
- **预处理指令**：引入必要的头文件。
- **宏定义**：定义矩阵大小\(n = 3\)。
```c
int main ()
{ int i,j,k,row,v,t;
int M[n];
float m,r,max0,max1;
float a[n][n],b[n],s[n],l[n][n],u[n][n],x[n],y[n];
a[0][0]=0.5,a[0][1]=1.1,a[0][2]=3.1;
a[1][0]=5.0,a[1][1]=0.96,a[1][2]=6.5;
a[2][0]=2.0,a[2][1]=4.5,a[2][2]=0.36;
b[0]=6.0,b[1]=0.96,b[2]=0.02;
```
- **变量声明**：声明了多个整型和浮点型数组用于存储矩阵和向量。
- **初始化**：初始化矩阵\(A\)和向量\(b\)。
##### 2. 分解阶段
```c
for(k=0;k<n;k++)
{ for(i=k;i<n;i++)
{ m=0;
if(k-1>=0)
{ for(t=0;t<k;t++)
m=m+l[i][t]*u[t][k];
}
s[i]=a[i][k]-m;
}
max0=abs(s[k]);
for(i=k;i<n;i++)
{ if(abs(s[i])>max0)
{ max0=abs(s[i]);
row=i;}
}
```
- **计算中间向量\(s\)**：\(s\)向量存储了未更新的主对角线元素。
- **选取主元**：找到绝对值最大的元素位置，并记录到`row`。
```c
max1=s[k];
M[k]=row;
if(row!=k)
{ if(k-1>=0)
{ for(t=0;t<=k-1;t++)
{ r=l[k][t];
l[k][t]=l[row][t];
l[row][t]=r;
}
}
for(t=k;t<n;t++)
{ r=a[k][t];
a[k][t]=a[row][t];
a[row][t]=r;
}
r=s[k];
s[k]=s[row];
s[row]=r;
}
```
- **交换行**：如果选取的主元不在当前列的顶部，则进行行交换操作。
- **更新\(s\)向量**：确保\(s\)向量始终表示主对角线元素。
```c
u[k][k]=s[k];
if(k<n-1)
{ for(j=k+1;j<n;j++)
{ m=0;
if(k-1>=0)
{ for(t=0;t<k;t++)
m=m+l[k][t]*u[t][j];
}
u[k][j]=a[k][j]-m;
}
for(i=k+1;i<n;i++)
l[i][k]=s[i]/u[k][k];
}
}
```
- **更新\(L\)和\(U\)矩阵**：根据选主元结果，逐步填充\(L\)和\(U\)矩阵。
##### 3. 求解线性方程组
```c
for(k=0;k<n-1;k++)
{ v=M[k];
m=b[k];
b[k]=b[v];
b[v]=m;
}
```
- **交换向量\(b\)中的元素**：根据主元选择时的行交换情况，相应地调整\(b\)向量。
```c
y[0]=b[0];
for(i=1;i<n;i++)
{ m=0;
for(t=0;t<=i-1;t++)
m=m+l[i][t]*y[t];
y[i]=b[i]-m;
}
```
- **前向代换**：利用\(L\)矩阵和\(b\)向量求得中间向量\(y\)。
```c
x[n-1]=y[n-1]/u[n-1][n-1];
for(i=n-2;i>=0;i--)
{ m=0;
for(t=i+1;t<n;t++)
m=m+u[i][t]*x[t];
x[i]=(y[i]-m)/u[i][i];
}
```
- **后向代换**：利用\(U\)矩阵和\(y\)向量求得最终解向量\(x\)。
```c
for(i=0;i<n;i++)
printf("%f\n",x[i]);
getch();
}
```
- **输出结果**：打印解向量\(x\)的每个元素。
#### 四、总结
通过上述解析可以看出，这段C语言程序实现了选主元的Doolittle分解法求解线性方程组的功能。该程序不仅考虑了数值稳定性的问题，还通过选主元策略提高了计算的精度。对于理解和掌握数值分析中的关键算法，尤其是涉及矩阵分解的方法来说，本程序提供了一个很好的示例。
qq_521503762022-04-07运行不了还在这发
