求C语言课程设计：用高斯列主元消元法解线性方程组

这里向你推荐一下克鲁特算法（其实就是对高斯列主元消元法进行优化，使之更适合于计算机编程），首先将矩阵A进行LU分解（将系数矩阵分解成一个上三角矩阵和一个下三角矩阵），分解的过程中用到了隐式的主元寻找法，同时利用克鲁特算法可以将两个nn矩阵压缩到一个nn矩阵中，大大节省了存储空间提高了计算速度。

方程可化为LUx=B，令Ux=y --->Ly=B

然后利用回代先求y，再利用y求x

因为该方法在求解过程中不涉及增广矩阵所以矩阵B几乎不参与什么运算，所以它的计算速度应该能够达到高斯列主元消元法的三倍，但原理与其基本一致。

而且我在程序中使用了动态数组方便你今后进行扩展。

以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序

如果你需要详细的理论讲解我可以将这两本书和源程序发给你，上面的论述相当详细足够你答辩用的了，我的邮箱hu_hu605@163com

计算结果：

A矩阵:

2 2 5

3 4 7

1 3 3

B矩阵:

5

6

5

解矩阵:

x 1=-7

x 2=0333333

x 3=366667

Press any key to continue

#include <cmath>

#include <iostream>

#include <iomanip>

#include <cstdlib>

#include <functional>

#include <vector>

#include <algorithm>

using namespace std;

#define TINY 10e-20 //A small number

#define N 3

void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d);//对矩阵进行LU分解

void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b);//对矩阵进行前向和后向回代

void root(vector<vector<float> > &x,vector<float> &col);//解方程结果保存在y中

void iniv(vector<vector<float> > &x,vector<float> line,int n);//对二维动态数组进行初始化

void main()

{

int i,j,n=N;//输入矩阵的维数

float A[N][N]={{2,2,5},{3,4,7},{1,3,3}};//左边A矩阵

float B[N]={5,6,5};//右边B矩阵

vector<vector<float> > x;//建立动态二维数组存放A，保证你的程序进行扩展时只改A，B，N

vector<float> line;

vector<float> y(n);//建立动态数组存放B

iniv(x,line,n);

yclear();

for(i=0;i<n;i++)//将A赋给x，B赋给y

{

ypush_back(B[i]);

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

{

x[i]push_back(A[i][j]);

}

}

cout<<"A矩阵:"<<endl;

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

{

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

{

cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";

}

cout<<endl;

}

cout<<"B矩阵:"<<endl;

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

{

cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;

}

root(x,y);//求根

cout<<"解矩阵:"<<endl;

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

{

cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl;

}

cout<<endl;

}

void root(vector<vector<float> > &x,vector<float> &col)

{

int n=xsize(),i=0,j=0;

vector<int> index(n);//用于记录寻找主元素过程中对矩阵的初等变换

indexclear();

float m=10;//记录变换方式，此程序中无用

ludcmp(x,n,index,m);//进行LU分解

lubksb(x,n,index,col);//根据分解结果进行回带

}

//以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序

//如果你需要详细的理论讲解我可以将这两本书和源程序发给你，我的邮箱hu_hu605@163com

void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)

{

int i,imax,j,k;

float big=0,dum=0,sum=0,temp=0;

vector<float> vv(n);

vvclear();

d=10;

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

{

big=00;

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

if ((temp=fabs(a[i][j])) > big)

big=temp;

vv[i]=10/big;

}

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

{

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

{

sum=a[i][j];

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

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

a[i][j]=sum;

}

big=00;

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

{

sum=a[i][j];

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

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

a[i][j]=sum;

if ( (dum=vv[i]fabs(sum)) >= big)

{

big=dum;

imax=i;

}

}

if (j != imax)

{

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

{

dum=a[imax][k];

a[imax][k]=a[j][k];

a[j][k]=dum;

}

d = -(d);

vv[imax]=vv[j];

}

indx[j]=imax;

if (a[j][j] == 00)

a[j][j]=TINY;

if (j != n)

{

dum=10/(a[j][j]);

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

a[i][j] = dum;

}

}

}

void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)

{

int i,ii=0,ip,j;

float sum;

for(i=0;i<n;i++)//按LU分解时寻找主元所进行的初等变换进行反边变换。

{

ip=indx[i];

sum=b[ip];

b[ip]=b[i];

b[i]=sum;

}

sum=0;

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

{

sum=0;

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

{

sum+=a[i][j]b[j];

}

b[i]=b[i]-sum;

}

b[n-1]=b[n-1]/a[n-1][n-1];

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

{

sum=0;

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

{

sum+=a[i][j]b[j];

}

b[i]=(b[i]-sum)/a[i][i];

}

}

void iniv(vector<vector<float> > &x,vector<float> line,int n)

{

int i,j;

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

{

xpush_back(line);

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

{

x[i]clear();

}

}

}

程序如下

function x=gauss(A,b) %高斯求解方程组

%x=gauss(A,b)

n=length(A);

a=[A,b];

for k=1:n-1

maxa=max(abs(a(k:n,k)));

if maxa==0

return;

end

for i=k:n

if abs(a(i,k))==maxa

y=a(i,k:n+1);a(i,k:n+1)=a(k,k:n+1);a(k,k:n+1)=y;

break;

end

end

for i=k+1:n

l(i,k)=a(i,k)/a(k,k);

a(i,k+1:n+1)=a(i,k+1:n+1)-l(i,k)a(k,k+1:n+1);

end

end

%回代

if a(n,n)==0

return

end

x(n)=a(n,n+1)/a(n,n);

for i=n-1:-1:1

x(i)=(a(i,n+1)-sum(a(i,i+1:n)x(i+1:n)))/a(i,i);

end

调用示例如下：

A=[2,-1,3;4,2,5;1,2,0];

b=[1;4;7];

x=gauss(A,b)

x =

9 -1 -6

列主元消去法

include "mathh"

#include "stdioh"

void zhuyuan (int k,int n,float a[20][21])

{int t,i,j;

float x,y;

x=fabs(a[k][k]);t=k;

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

if (fabs(a[i][k])>x)

{x=fabs(a[i][k]);t=i;}

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

{y=a[k][j];a[k][j]=a[t][j];a[t][j]=y;}

}

void xiaoyuan(int n,float a[20][21])

{int k,i,j;

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

{zhuyuan(i,n,a);

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

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

a[j][k]=a[j][k]-a[j][i]a[i][k]/a[i][i];

}

}

void huidai(int n,float a[20][21],float x[20])

{int i,j;

x[n]=a[n][n+1]/a[n][n];

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

{ x[i]=a[i][n+1];

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

x[i]=x[i]-a[i][j]x[j];

x[i]=x[i]/a[i][i];

}

}

void main ()

{int n,i,j;float a[20][21],m[20];

printf("please input n=");

scanf("%d",&n);

n=n-1;

printf("please input aij and bi:\n");

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

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

{if(j==n+1)

printf("b%d=",i+1);

else

printf("a%d%d=",i+1,j+1);

scanf("%f",&a[i][j]);

}printf("\n");

}

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

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

if(j==n+1)

printf("%f\n",a[i][j]);

else

printf("%f ",a[i][j]);

}

xiaoyuan(n,a);

huidai(n,a,m);

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

printf("x%d=%f\n",i+1,m[i]);

printf("\n");

}

最小二乘法

#include "stdioh"

#include "conioh"

#include "stdlibh"

#include "mathh"

#define N 9//N个节点

#define M 2//M次拟合

#define K 2M

void zhuyuan (int k,int n,float a[M+1][M+2])

{int t,i,j;

float x,y;

x=fabs(a[k][k]);t=k;

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

if (fabs(a[i][k])>x)

{x=fabs(a[i][k]);t=i;}

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

{y=a[k][j];a[k][j]=a[t][j];a[t][j]=y;}

}

void xiaoyuan(int n,float a[M+1][M+2])

{int k,i,j;

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

{zhuyuan(i,n,a);

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

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

a[j][k]=a[j][k]-a[j][i]a[i][k]/a[i][i];

}

}

void huidai(int n,float a[M+1][M+2],float x[M+1])

{int i,j;

x[n]=a[n][n+1]/a[n][n];

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

{ x[i]=a[i][n+1];

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

x[i]=x[i]-a[i][j]x[j];

x[i]=x[i]/a[i][i];

}

}

void main()

{float x_y[N][2],A[N][K+1],B[N][M+1],AA[K+1],BB[M+1],a[M+1][M+2],m[M+1];

int i,j,n;

printf("请输入%d个已知点：\n",N);

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

{

printf("(x%d y%d):",i,i);

scanf("%f %f",&x_y[i][0],&x_y[i][1]);

}

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

{

A[i][0]=1;

for(j=1;j<=K;j++)

A[i][j]=A[i][j-1]x_y[i][0];

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

B[i][j]=A[i][j]x_y[i][1];

}

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

for(AA[j]=0,i=0;i<N;i++)

AA[j]+=A[i][j];

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

for(BB[j]=0,i=0;i<N;i++)

BB[j]+=B[i][j];

for(i=0;i<M+1;i++)

{a[i][M+1]=BB[i];

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

a[i][j]=AA[i+j];}

n=M;

printf("正规系数矩阵为:\n");

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

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

printf("%f ",a[i][j]);

printf("\n");}

xiaoyuan(n,a);

huidai(n,a,m);

printf("拟合曲线方程为:\ny(x)=%g",m[0]);

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

{

printf(" + %g",m[i]);

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

{

printf("X");

}

}

}

库塔

#include "stdioh"

float f(float x,float y)

{float p;

p=x+y;

return p;

}

void main()

{float x,y,k1,k2,k3,k4,h,a,b;

int i,n;

printf("请输入区间:");

scanf("%f,%f",&a,&b);

printf("(%f,%f)\n",a,b);

printf("请输入步长h=");

scanf("%f",&h);

printf("请输入y的初值:");

scanf("%f",&y);

x=a;

printf("x=%f,y=%f\n",x,y);

n=(b-a)/h;

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

{k1=hf(x,y);

k2=hf(x+h/2,y+k1/2);

k3=hf(x+h/2,y+k2/2);

k4=hf(x+h,y+k3);

y+=(k1+2k2+2k3+k4)/6;

x=a+hi;

printf("x=%f,y=%f\n",x,y);

}

}

贝格

#include <stdioh>

#include <mathh>

#include <conioh>

double f(double x)

{double p;

p=1/x;

return p;}

void main()

{double T[6],S[5],C[4],R[3],ff,b,a,e;

int i,j,k;

scanf("%lf,%lf,%lf",&a,&b,&e);

T[0]=(b-a)1/2(f(a)+f(b));

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

{ff=0;

for(i=1;2i<=pow(2,k);i++)

ff+=f(a+(b-a)/pow(2,k)(2i-1));

T[k]=1/2T[k-1]+(b-a)/pow(2,k-1)ff;}

for(i=0,j=-1,k=-2;i<5;i++,j++,k++)

{S[i]=(4T[i+1]-T[i])/3;

if(i>0)

C[j]=(16S[i]-S[i-1])/15;

if(j>0)

R[k]=(64C[j]-C[j-1])/63;

if(k>0)

if((R[k]-R[k-1])<e)

break;}

printf("%lf",R[k]);

}

以上就是关于求C语言课程设计：用高斯列主元消元法解线性方程组全部的内容，包括:求C语言课程设计：用高斯列主元消元法解线性方程组、matlab 用列主元高斯消去法求逆矩阵、列主元消去法 最小二乘法 龙贝格公式 标准四阶龙格库塔的编程等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！