高斯消元初步(Gauss算法)
Gauss算法,称为高斯消元算法,用来解决n元一次方程,在解决线性方程问题起着重要作用。
简述
运用高斯消元的方法,我们可以在O(n3)的时间求出n元线性方程,但是由于时间复杂度的原因,请注意题目数据范围的提示。
高斯消元三大定理(在小学就学过了吧):
1.两个方程互换位置,解不变;
2.一个方程进行加减乘除,解不变;
3.一个方程乘上数k加上另一个方程,解不变;
这便是我们解决的基础;
过程:
这里给出luogu例题链接,这样方便寻找;
我们这里不用luogu的样例示范(因为不是整数好麻烦),这里给出方程
2 x + 3 y - z = 21;
x + 2 y + 2 z = 7;
3 x + y + 5 z = 8;
我们将系数提出,然后就可以得到一个3 * 3的矩阵,之后将每个方程等号右边放到矩阵的最右边,就得到了:
这里每个方程的结果与系数我用黑线隔开了,想必也更清楚;
有了定理,我们理一下目标:
我们如果将每一个方程只留下一个未知数的系数,那么最后就可以求解了,如:
当然系数不一定只会是1,但是只要除一下就好,根据这个定义,我们将第 i 个未知数的前系数非零而且其他系数都为零,这个系数在矩阵的位置为 i ,i;
这样的矩阵称为“简化阶梯矩阵”;
我们只要将每个矩阵化成简化阶梯矩阵即可;
步骤:
1.枚举第 i 个未知数(外循坏);
2.决定在哪一行求解这个未知数:
这里采用先对每一行 j 第 i 个系数找到最大值,有最大值的这一行定义为第ms行(名字随便起的,没有其他意思),然后将第ms行交换至第 i 行
3.判断第i行第i个数的值是否为0,这里由于数学期望和精度问题,我们将这个判断改为这个数的值是否小于我们定义的那个精度,如果小于(那就相当于为0了),
那么无解(因为这个项的系数是所有中最大的,所以其他的也都为0,一定无解);
4.然后进行消元,就是将其他方程这个项的系数归0,这里有精度问题,但是从期望来讲,是不成问题的;
Code:
(我才不会说其实我有模拟操作但是太麻烦不想写了。。。)
不过我很良心,所以我有输出模拟,运行一下我的代码就行了;
#include<bits/stdc++.h>
#define maxn 107
#define db double
using namespace std;
int n;
const db cmp=1e-8;
db a[maxn][maxn]; //模拟啦 biu~
void biu(int x){
printf("work %d \n",x);
for(int i=1;i<=n;i++){
for(int j=1;j<=n+1;j++)
cout<<a[i][j]<<" ";
cout<<endl<<endl;
}
} int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
scanf("%lf",&a[i][j]);
for(int i=1;i<=n;i++){
int ms=i;
for(int j=i;j<=n;j++)
if(fabs(a[j][i])>fabs(a[ms][i]))
ms=j;
if(fabs(a[ms][i])<cmp){
puts("No Solution");
return 0;
}
for(int j=1;j<=n+1;j++)
swap(a[i][j],a[ms][j]);
biu(2);
for(int j=1;j<=n;j++){
if(j==i) continue;
db rate=a[j][i]/a[i][i];
for(int k=i+1;k<=n+1;k++)
a[j][k]-=a[i][k]*rate;
}
biu(4);//良心模拟。。。
}
for(int i=1;i<=n;i++) printf("%.2f\n",a[i][n+1]/a[i][i]);
return 0;
}
例题
题目描述
有一个球形空间产生器能够在 nnn 维空间中产生一个坚硬的球体。现在,你被困在了这个 nnn 维球体中,你只知道球面上 n+1n+1n+1 个点的坐标,你需要以最快的速度确定这个 nnn 维球体的球心坐标,以便于摧毁这个球形空间产生器。
输入格式
第一行是一个整数 nnn (1<=N=10)(1<=N=10)(1<=N=10)。接下来的 n+1n+1n+1 行,每行有 nnn 个实数,表示球面上一点的 nnn 维坐标。每一个实数精确到小数点后 666 位,且其绝对值都不超过 200002000020000。
输出格式
有且只有一行,依次给出球心的 nnn 维坐标( nnn 个实数),两个实数之间用一个空格隔开。每个实数精确到小数点后 333 位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。
这个只要构造出矩阵即可;
可惜我不会在博客上用数学公式。。。
所以就随便写写了;
设xi为第i维的坐标;
sum(j=1,n) { (a[ i , j ] - x [ j ]) }=r^2;
这样的方程共有11个,我们要将r消掉,所以将相邻的两个方程相减,得到11个方程,然后将多项式拆开,合并,移项得到;
sum(j=1,n){ 2*(a[ i , j ] - a[ i + 1 , j ) * x [ j ] } = sum(j=1,n){ a[ i , j ] - a[ i + 1 , j ] };
这样就可以将左边的作为方程左边,右边作为结果,列出矩阵,这里还不需要检验,直接上代码。。。
#include<bits/stdc++.h>
#define db double
using namespace std;
int n;
db a[17][17],c[17][17];
const db cmp=1e-8;
db calc(db x){return x*x;} int main(){
scanf("%d",&n);
for(int i=1;i<=n+1;i++){
for(int j=1;j<=n;j++){
scanf("%lf",&a[i][j]);
if(i==1) continue;
c[i-1][j]=2*(a[i-1][j]-a[i][j]);
c[i-1][n+1]+=calc(a[i-1][j])-calc(a[i][j]);
}
}
for(int i=1;i<=n;i++){
int ms=i;
for(int j=i+1;j<=n;j++)
if(fabs(c[j][i])>fabs(c[ms][i])) ms=j;
if(ms!=i) for(int j=1;j<=n+1;j++)
swap(c[i][j],c[ms][j]);
for(int j=1;j<=n;j++){
if(i==j) continue;
db rate=c[j][i]/c[i][i];
for(int k=i+1;k<=n+1;k++)
c[j][k]-=c[i][k]*rate;
}
}
for(int i=1;i<=n;i++) printf("%.3f ",c[i][n+1]/c[i][i]);
return 0; }
其实还有一些拓展内容,到时候再补充。。。
最新文章
- [Java 安全]消息摘要与数字签名
- 【Visual Lisp】人机交互与数据处理(表除外)-lisp
- C++问题-Qt Visual Studio Add-in
- Android测试分析3
- linux mysql 授权以及 iptables开启3306
- Android的minSdkVersion,targetSdkVersion,maxSdkVersion
- 简单实现div+css页面自适应
- 基于BrokerPattern服务器框架
- VBS常用使用技巧
- 每天一道Java题[7]
- 转:【Java并发编程】之二十三:并发新特性—信号量Semaphore(含代码)
- Function Programming - 纯函数(Pure Function)
- Linux下批量修改文件名方法
- Oracle解决Ora-01653无法扩展表空间问题
- redux源码解析-redux的架构
- EF 中多个数据库迁移
- Exception 01 : org.hibernate.engine.jndi.JndiException: Error parsing JNDI name [foo]
- HDU4283(KB22-G)
- 第五周可执行代码 以及 Canvas 制作个人PSP分类饼图
- Ios8 Xcode6 设置Launch Image 启动图片
热门文章
- rm: cannot remove `/tmp/localhost-mysql_cacti_stats.txt&#39;: Operation not permitted
- 【ORA】ORA-27101快速处理方法
- 两行代码修复了解析MySQL8.x binlog错位的问题!!
- dubbo快速入门demo
- 24V降压5V芯片,5A,4.5V-30V输入,同步降压调节器
- 低功耗降线性稳压器,24V转5V降压芯片
- 前端面试准备笔记之JavaScript(02)
- 可视化Go内存管理
- HDFS HBase Solr Which one? 从访问模式角度决策
- autoreload 线程 进程管理 并发的处理方法