说明：这是数值分析前几章算法的C++实现，可以解多元一次方程，用惯了matlab，有兴趣试试C++实现么？这是我做的，欢迎找BUG。我测试得的解均正确。不做代码说明了（函数划分写得非常清楚），具体可参考你的数值分析课本，再映照我写的函数。
// linear_equation.cpp: 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include<iostream>
#include<cmath>
#include<utility>
#include<string>
#include <sstream>
using namespace std;
void num_6_transefer(double num);
int zhanzhuan(int x, int y)
{
int max, min;
if(x>y)
{
max = x;
min = y;
}
else
{
max = y;
min = x;
}
int temp = max % min;
while (temp != 0)
{
max = min;
min = temp;
temp = max % min;
}
return min;
}
void Doolittle(int n)
{
double **p, sum;
int i, j, k;
p = new double *[n + 1];
for (i = 1; i <= n; ++i)
p[i] = new double[n + 2];
cout << "请输入方程的增广矩阵：" << '\n';
for (i = 1; i <= n; ++i)
for (j = 1; j <= n + 1; ++j)
cin >> p[i][j];
for (i = 2; i <= n; ++i)
{
p[i][1] /= p[1][1];
}
for (i = 2; i <= n; ++i)
{
for (j = i; j <= n; ++j)
{
sum = 0;
for (k = 1; k < i; ++k)
{
sum += p[i][k] * p[k][j];
}
p[i][j] -= sum;
}
for (j = i + 1; j <= n; ++j)
{
sum = 0;
for (k = 1; k < i; ++k)
{
sum += p[j][k] * p[k][i];
}
p[j][i] -= sum;
p[j][i] /= p[i][i];
}
}
cout << "该矩阵分解的下三角矩阵 L为：" << '\n';
cout << '1' << '\n';
for (i = 2; i <= n; ++i)
{
for (j = 1; j < i; ++j)
{
num_6_transefer(p[i][j]);
cout << ' ';
}
cout << '1' << '\n';
}
cout << "该矩阵分解的上三角矩阵U为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j<i; ++j)
cout << ' ';
for (j = i; j <= n; ++j)
{
num_6_transefer(p[i][j]);
cout << ' ';
}
cout << '\n';
}
for (i = 1; i <= n; ++i)
{
for (j = 1; j < i; ++j)
{
p[i][n + 1] -= p[i][j] * p[j][n + 1];
}
}
for (i = n; i >= 1; --i)
{
for (j = i + 1; j <= n; ++j)
{
p[i][n + 1] -= p[i][j] * p[j][n + 1];
}
p[i][n + 1] /= p[i][i];
}
cout << "用Doolittle法的解为：" << '\n';
for (i = 1; i <= n; ++i)
{
num_6_transefer(p[i][n+1]);
cout << endl;
}
for (i = 1; i <= n; ++i)
delete[]p[i];
delete[]p;
}
void num_6_transefer(double num)//对解出的解进行处理精确到小数点后六位,化为分数同时对分数进行约分,并约分
{
pair<int, int> fraction(-101, -101);
if (abs(num - 0) < 0.000001)
{
fraction.first = 0;
}
else
{
double i, j;
for (i = -100; i <= 100; i += 1)
for (j = -100; j < 0; ++j)
{
if (abs(num - i / j) < 0.000001)
{
fraction.first = (int)i / zhanzhuan(i, j);
fraction.second = (int)j / zhanzhuan(i, j);
}
}
for (i = -100; i <= 100; i += 1)
for (j = 1; j < 100; ++j)
{
if (abs(num - i / j) < 0.000001)
{
fraction.first = (int)i / zhanzhuan(i, j);
fraction.second = (int)j / zhanzhuan(i, j);
}
}
}
if (fraction.first == -101)
cout << num;
else if (fraction.first == 0)
cout << '0' ;
else if (fraction.second == 1)
cout << fraction.first;
else if (fraction.second == -1)
cout <<'-'<< fraction.first;
else if (fraction.first != fraction.second)
{
if (fraction.second < 0)
{
fraction.first *= -1;
fraction.second *= -1;
}
cout << fraction.first << "/" << fraction.second ;
}
else
cout << '1' ;
}
void GS_iteration(int n)
{
double **p,*x,*temp, sum;
int i, j,k;
bool flag;
p = new double *[n + 1];
x = new double[n + 1];
temp = new double[n + 1];
for (i = 1; i <= n; ++i)
p[i] = new double[n + 2];
cout << "请输入方程的增广矩阵：" << '\n';
for (i = 1; i <= n; ++i)
for (j = 1; j <= n + 1; ++j)
cin >> p[i][j];
for (i = 1; i <= n; ++i)
{
temp[i]=x[i] = 0;
}
for(i=1;1;++i)
{flag = 1;
for (j = 1; j <= n; ++j)
{
sum = 0;
for (k = 1; k <j; ++k)
sum += p[j][k] * x[k];
for (k = j+1; k <=n; ++k)
sum += p[j][k] * x[k];
x[j] = 1 / p[j][j] * (p[j][n + 1] - sum);
}
for (j = 1; j <= n; ++j)
{
if ((x[j] - temp[j])>=0.00001|| (temp[j] - x[j]) >= 0.00001)
{
temp[j] = x[j];
//cout << "xxx xxx ";
flag = 0;
}
}
if (flag)
{
cout << "最终迭代结果为：" << '\n';
for (j = 1; j <= n; ++j)
{
cout << 'x' << j << '=' << x[j] << ' ';
}
cout << '\n';
return;
}
cout << "第"<<i<<"次迭代结果为：" << '\n';
for (j = 1; j <= n; ++j)
{
cout << 'x' << j << '=' << x[j] << ' ';
}
cout<< '\n';
}
for (i = 1; i <= n; ++i)
delete[]p[i];
delete[]x;
delete[]temp;
}
double** matrix_inversion(double **p, int n )
{
double **result,**ni,temp;
int i, j, k;
result = new double *[n + 1];
for (i = 1; i <= n; ++i)
result[i] = new double[2*n + 1];
ni = new double *[n + 1];
for (i = 1; i <= n; ++i)
ni[i] = new double[ n + 1];
for (i = 1; i <= n; ++i)
for (j = 1; j <=2 *n; ++j)
{ if(j<=n)
result[i][j] = p[i][j];
else if(j!=n+i)
result[i][j] = 0;
else
result[i][j] = 1;
}
for (j = 1; j < n; ++j)//j代表列
{
for (i = j; i <= n; ++i)
{
if (result[i][j] != 0)
break;
}
if (i != j)
{
for (k = j; k <= 2 * n; ++k)
swap(result[i][k], result[j][k]);
}
for (i = j + 1; i <= n; ++i)
{
if (result[i][j] != 0)
{
temp = result[i][j] / result[j][j];
result[i][j] = 0;
for (k = j + 1; k <= 2 * n; ++k)
result[i][k] -= result[j][k] * temp;
}
}
}
for (i = n+1; i <= 2 * n; ++i)
result[n][i] /= result[n][n];
result[n][n] = 1;//将最后一行置为1
for (j = n; j>1; --j)
{
for (i = j - 1; i >= 1; --i)
{
if (result[i][j] != 0)
{
for (k = j + 1; k <= 2 * n; ++k)
result[i][k] -= result[j][k]*result[i][j] ;
result[i][j] = 0;
}
}
for (k = j + 1; k <= 2 * n; ++k)
result[j - 1][k] /= result[j - 1][j - 1];
result[j - 1][j - 1] = 1;//为上一列置1
}
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= n; ++j)
ni[i][j] = result[i][n + j];
}
for (i = 1; i <= n; ++i)
delete[]result[i];
delete[]result;
return ni;
}
double ** matrix_mul(double **p, double **q,int n)
{
double **result;
int i, j,k;
result = new double *[n + 1];
for (i = 1; i <= n; ++i)
result[i] = new double[n + 1];
for(i=1;i<=n;++i)
for (j = 1; j <= n; ++j)
{
result[i][j] = 0;
for (k = 1; k <= n; ++k)
{
result[i][j]+= p[i][k] * q[k][j];
}
}
for (i = 1; i <= n; ++i)
delete []p[i];
delete[]p;
return result;
}
void matrix_ana(int n)
{
double **p, **pt,**pni, sum, max1 = 0,maxpni1=0, maxmax= 0, maxF=0;
int i, j;
p = new double *[n + 1];
for (i = 1; i <= n; ++i)
p[i] = new double[n + 1];
pt = new double *[n + 1];
for (i = 1; i <= n; ++i)
pt[i] = new double[n + 1];
cout << "请输入方程的系数矩阵：" << '\n';
for (i = 1; i <= n; ++i)
{
sum = 0;
for (j = 1; j <= n; ++j)
{
cin >> p[i][j];
maxF += p[i][j] * p[i][j];
if (p[i][j] > 0)
sum += p[i][j];
else
sum -= p[i][j];
}
if(sum>maxmax)
maxmax =sum;
}
maxF = sqrt(maxF);
for (i = 1; i <= n; ++i)
{
sum = 0;
for (j = 1; j <= n; ++j)
{
if (p[j][i] > 0)
sum += p[j][i];
else
sum -= p[j][i];
}
if (sum>max1)
max1 = sum;
}
for(i=1;i<=n;++i)
for (j = 1; j <= n; ++j)
{
pt[j][i] = p[i][j];
}
cout << "该矩阵的转置矩阵为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= n; ++j)
cout << pt[i][j] << ' ';
cout << '\n';
}
pni = matrix_inversion(p, n);
cout << "该矩阵的逆矩阵为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= n; ++j)
{
num_6_transefer(pni[i][j]);
cout << ' ';
}
cout << '\n';
}
cout << "该矩阵的1-范数为：" <<max1<< '\n';
cout << "该矩阵的∞-范数为：" << maxmax << '\n';
cout << "该矩阵的F-范数为：" << maxF << '\n';
pt = matrix_mul(pt, p, n);
cout << "该矩阵的ATA为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= n; ++j)
cout << pt[i][j]<<' ';
cout << '\n';
}
for (i = 1; i <= n; ++i)
{
sum = 0;
for (j = 1; j <= n; ++j)
{
if (pni[j][i] > 0)
sum += pni[j][i];
else
sum -= pni[j][i];
}
if (sum>maxpni1)
maxpni1 = sum;
}
cout << "该矩阵的1条件数为："<<max1*maxpni1 << '\n';
for (i = 1; i <= n; ++i)
{
delete[]p[i];
}
delete[]p;
delete[]pni;
}
void Crout(int n)
{
n = 3;
int i, j, button;
double *a, *b, *c, *d;
a = new double[n + 1];
b = new double[n+1];
c = new double[n];
d = new double[n+1];
cout << n;
cout << "请输入方程的增广矩阵：" << '\n';
for (i = 1; i <= n; ++i)
for (j = 1; j <= n + 1; ++j)
{
if (j == i)
cin >> a[i];
else if ((j == i + 1)&&(i!=n))
cin >> c[i];
else if (j == i - 1)
cin >> d[i];
else if (j == n + 1)
cin >> b[i];
else
cin >> button;
}
c[1]/= a[1];
for (i = 2; i <= n; ++i)
{
a[i] -= d[i] * c[i - 1];
c[i] /= a[i];
}
cout << "所分解矩阵T为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= n; ++j)
{
if (j == i)
{
num_6_transefer(a[i]);
cout << ' ';
}
else if (j == i - 1)
{
num_6_transefer(d[i]);
cout << ' ';
}
else
cout << '0' << ' ';
}
cout << '\n';
}
cout << "所分解矩阵M为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= n; ++j)
{
if (j == i)
cout << 1 << ' ';
else if (j == i + 1)
cout << c[i] << ' ';
else
cout << '0' << ' ';
}
cout << '\n';
}
b[1] /= a[1];
for (i = 2; i <= n; ++i)
{
b[i] -= b[i - 1] * d[i];
b[i] /= a[i];
}
for (i = n - 1; i >= 1; --i)
b[i] -= c[i] * b[i + 1];
cout << "用Crout法的解为：" << '\n';
for (i = 1; i <= n; ++i)
{
cout << b[i] << '\n';
}
delete[]a;
delete[]b;
delete[]c;
delete[]d;
}
void Cholesky(int n)
{
double **p, button, sum = 0;
int i, j, k;
p = new double *[n + 1];
for (i = 1; i <= n; ++i)
p[i] = new double[i + 1];
cout << "请输入方程的增广矩阵：" << '\n';
for (i = 1; i <= n; ++i)
for (j = 1; j <= n + 1; ++j)
{
if (j <= i)
cin >> p[i][j];
else if (j < n + 1)
cin >> button;
else
cin >> p[i][i + 1];
}
p[1][1] = sqrt(p[1][1]);
for (i = 2; i <= n; ++i)
{
p[i][1] /= p[1][1];
}
for (i = 2; i <= n; ++i)
{
sum = 0;
for (j = 1; j < i; ++j)
{
sum += p[i][j] * p[i][j];
}
p[i][i] -= sum;
p[i][i] = sqrt(p[i][i]);
for (j = i + 1; j <= n; ++j)
{
sum = 0;
for (k = 1; k < i; ++k)
{
sum += p[i][k] * p[j][k];
}
p[j][i] -= sum;
p[j][i] /= p[i][i];
}
}
cout << "按平方根法分解所得矩阵G为：" << '\n';
for (i = 1; i <= n; ++i)
{
for (j = 1; j <= i; ++j)
{
if (p[i][j]>0)
cout << "√";
else
cout <<'-'<< "√";
num_6_transefer(pow(p[i][j], 2));
cout << ' ';
}
cout << '\n';
}
p[1][2] /= p[1][1];
for (i = 2; i <= n; ++i)
{
for (j = 1; j < i; ++j)
{
p[i][i + 1] -= p[i][j] * p[j][j + 1];
}
p[i][i + 1] /= p[i][i];
}
p[n][n + 1] /= p[n][n];
for (i = n - 1; i >= 1; --i)
{
for (j = n; j > i; --j)
{
p[i][i + 1] -= p[j][i] * p[j][j + 1];
}
p[i][i + 1] /= p[i][i];
}
cout << "用Cholesky法的解为：" << '\n';
for (i = 1; i <= n; ++i)
{
num_6_transefer(p[i][i + 1]);
cout << endl;
}
for (i = 1; i <= n; ++i)
delete []p[i];
delete[]p;
}
void list_gauss(int n)
{
double **p, l,max,now;
int i, j, k,maxpos;
p = new double *[n + 1];
for (i = 1; i <= n; ++i)
p[i] = new double[n + 2];
cout << "请输入方程的增广矩阵：" << '\n';
for (i = 1; i <= n; ++i)
for (j = 1; j <= n + 1; ++j)
cin >> p[i][j];
for (i = 1; i <= n - 1; ++i)
{
now=max = p[i][i];
maxpos = 1;
if (max < 0)
max *= -1;
for (j = i + 1; j <= n; ++j)
{
if (p[j][i] < 0)
now = p[j][i] * -1;
else
now = p[j][i];
if (now > max)
{
max = now;
maxpos = j;
}
}
if (max != p[i][i])
{
for (k = i; k <= n + 1; ++k)
{
swap(p[maxpos][k], p[i][k]);
}
}
for (j = i + 1; j <= n; ++j)
{
l = p[j][i] / p[i][i];
for (k = i; k <= n + 1; ++k)
{
p[j][k] -= p[i][k] * l;
}
}
}
cout << "系数矩阵消元后为：";
for (i = 1; i <= n; ++i)
{
cout << '\n';
for (j = 1; j <= n + 1; ++j)
{
num_6_transefer(p[i][j]);
cout << ' ';
}
}
cout << endl;
for (i = n; i >= 1; --i)
{
for (j = i+1; j <=n; ++j)
{
p[i][n + 1] -= p[i][j] * p[j][n + 1];
}
p[i][n + 1] /= p[i][i];
}
cout<< "用列主元消去法的解为："<<'\n';
for (i = 1; i <= n; ++i)
{
num_6_transefer(p[i][n + 1]);
cout << endl;
}
for (i = 1; i <= n; ++i)
delete[]p[i];
delete[]p;
}
void gauss(int n)
{
double **p, l;
int i, j, k;
p = new double *[n + 1];
for (i = 1; i <= n; ++i)
p[i] = new double[n + 2];
cout << "请输入方程的增广矩阵：" << '\n';
for (i = 1; i <= n; ++i)
for (j = 1; j <= n + 1; ++j)
cin >> p[i][j];
for (i = 1; i <= n - 1; ++i)
for (j = i + 1; j <= n; ++j)
{
l = p[j][i] / p[i][i];
for (k = i; k <= n + 1; ++k)
{
p[j][k] -= p[i][k] * l;
}
}
cout << "系数矩阵消元后为：";
for (i = 1; i <= n; ++i)
{
cout << '\n';
for (j = 1; j <= n + 1; ++j)
{
num_6_transefer(p[i][j]);
cout << ' ';
}
}
cout << endl;
for (i = n; i >= 1; --i)
{
for (j = i + 1; j <= n; ++j)
{
p[i][n + 1] -= p[i][j] * p[j][n + 1];
}
p[i][n + 1] /= p[i][i];
}
cout << '\n' << "用顺序高斯法方程的解为：" << '\n';
for (i = 1; i <= n; ++i)
{
num_6_transefer(p[i][n + 1]);
cout << endl;
}
for (i = 1; i <= n; ++i)
delete[]p[i];
delete[]p;
}
int main()
{
char n;
while (1)
{
system("cls");
cout << "1:顺序高斯消去法:\n";
cout << "2:行列高斯消去法:\n";
cout << "3:直接三角分解法:\n";
cout << "4:正定方程的平方根算法:\n";
cout << "5:追赶法求解三对角方程:\n";
cout << "6:对方程组形态分析:\n";
cout << "7:使用高斯迭代法对方程组近似求解:\n";
cout << "8:表达式求值:\n";
cin >> n;
switch (n)
{
case '1':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
gauss(n);
break;
case '2':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
list_gauss(n);
break;
case '3':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
Doolittle(n);
break;
case '4':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
Cholesky(n);
break;
case '5':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
Crout(n);
break;
case '6':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
matrix_ana(n);
break;
case '7':
system("cls");
cout << "请输入方程的阶数：";
cin >> n;
GS_iteration(n);
break;
default:
cout << "输入有误";
}
cout << "按任意键继续";
cin >> n;
}
return 0;
}
298
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
