哪位高手把解线形方程的数值算法程序贴一下，小弟拜谢

shanglei · 发表于 2004-8-3 17:38:04

因为要做实习，哪位高手把解线形方程的数值算法程序贴一下，小弟拜谢

liang261898 · 发表于 2004-8-4 00:39:04

<

> 求解形如AX=B的线形方程组有直接法和迭代法，直接法有高斯消元法和列主元消元法。迭代法有雅可比迭代法和高斯塞德尔迭代法。<

> 雅可比迭代法程序高斯塞德尔迭代法程序<

>function s=jacobi(a,b,x0,eps) function s=gauss(a,b,x0,eps) D=diag(diag(a)); D=diag(diag(a)); D=inv(D); L=tril(a,-1); L=tril(a,-1); u=triv(b,1); u=triv(a,1); B=inv(D+L); B=-D*(L+u); f=c*b f=D*b; s=B*x0+f; s=B*x0+f; while norm(s-x0)>=eps while norm(s-x0)>=eps x0=s; x0=s; s=B*x0+f; s=B*x0+f; end end return return

knight · 发表于 2004-8-4 09:40:34

<

>//Jacobi 方法,以下是完整的c++程序，把矩阵A存在同目录下的A.txt中，向量b存在同目录下的B.txt中，结果在out.txt中，解得是Ax=b<

>#include <vector>
#include <math.h>
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h><

>using namespace std;ofstream out("out.txt");double operator * (vector<double> a,vector<double> b)
{
 double data=0;
 if(a.size()==b.size())
 {
 vector<double> :: iterator it1=a.begin();
 vector<double> :: iterator it2=b.begin();
 for(;it1!=a.end();++it1,++it2)
 {
 data+=((*it2)*(*it1));
 }
 }
 return data;
}double Delta(vector<double> a,vector<double> b)
{
 double data=0;
 if(a.size()==b.size())
 {
 vector<double> :: iterator it1=a.begin();
 vector<double> :: iterator it2=b.begin();
 for(;(it1!=a.end())&&(it2!=b.end());++it1,++it2)
 {
 if(fabs((*it1)-(*it2))>data)
 {
 data=fabs((*it1)-(*it2));
 }
 }
 }
 return data;
}

vector<double> Jacobi(vector<vector<double> > A,vector<double> B,double delta)
{
 vector<double> x1(B.size(),0);
 vector<double> x2(B.size(),0);
 do
 {
 x2=x1;
 for(int i=0;i<x2.size();++i)
 {
 double xx=A*x2;
 xx-=A*x2;
 x1=(B-xx)/A;
 }
 }while(Delta(x1,x2)>delta);
 return x1;
}void Print(vector<double> a)
{
 vector<double> :: iterator it=a.begin();
 for(;it!=a.end();++it)
 {
 cout<<(*it)<<"\t";
 out<<(*it)<<"\t";
 }
 cout<<endl;
 out<<endl;
}void Print(vector<vector<double> > A)
{
 for(int i=0;i<A.size();++i)
 {
 Print(A);
 }
}

int main()
{
 ifstream inA("A.txt");
 ifstream inB("B.txt");
 cout<<"Define the dimension of the matrix:"<<endl;
 int num;
 cin>>num;
 vector<vector<double> > A(num);
 for(int i=0;i<num;++i)
 {
 double b=0;
 for(int j=0;j<num;++j)
 {
 {
 A.push_back(1/((double)i+(double)j+1));
 b+=(1/((double)i+(double)j+1));
 }
 }
 B.push_back(b);
 }
 cout<<"Input the delta:"<<endl;
 double delta;
 cin>>delta;
 cout<<"The Result is:"<<endl;
 Print(Jacobi(A,B,delta));
 system("PAUSE");
}

knight · 发表于 2004-8-4 09:43:25

<

>//高斯塞德尔迭代法，用法同上<

>#include <vector>
#include <math.h>
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h><

>using namespace std;ofstream out("out.txt");double operator * (vector<double> a,vector<double> b)
{
 double data=0;
 if(a.size()==b.size())
 {
 vector<double> :: iterator it1=a.begin();
 vector<double> :: iterator it2=b.begin();
 for(;it1!=a.end();++it1,++it2)
 {
 data+=((*it2)*(*it1));
 }
 }
 return data;
}double Delta(vector<double> a,vector<double> b)
{
 double data=0;
 if(a.size()==b.size())
 {
 vector<double> :: iterator it1=a.begin();
 vector<double> :: iterator it2=b.begin();
 for(;(it1!=a.end())&&(it2!=b.end());++it1,++it2)
 {
 if(fabs((*it1)-(*it2))>data)
 {
 data=fabs((*it1)-(*it2));
 }
 }
 }
 return data;
}

vector<double> GS(vector<vector<double> > A,vector<double> B,vector<vector<double> > D,vector<vector<double> > L,double delta)
{
 vector<double> x1(B.size(),0);
 vector<double> x2(B.size(),0);
 do
 {
 x1=x2;
 for(int i=0;i<x2.size();++i)
 {
 double xx=L*x2;
 xx+=D*x1;
 x2=(B-xx)/A;
 }
 }while(Delta(x1,x2)>delta);
 return x1;
}void Print(vector<double> a)
{
 vector<double> :: iterator it=a.begin();
 for(;it!=a.end();++it)
 {
 cout<<(*it)<<"\t";
 out<<(*it)<<"\t";
 }
 cout<<endl;
 out<<endl;
}void Print(vector<vector<double> > A)
{
 for(int i=0;i<A.size();++i)
 {
 Print(A);
 }
}

int main()
{
 cout<<"Define the dimension of the matrix:"<<endl;
 int num;
 cin>>num;
 vector<vector<double> > A(num);
 vector<double> B;
 for(int i=0;i<num;++i)
 {
 double bb=0;
 for(int j=0;j<num;++j)
 {
 {
 A.push_back(1/((double)i+(double)j+1));
 bb+=(1/((double)i+(double)j+1));
 }
 }
 B.push_back(bb);
 }
 cout<<"Input the delta:"<<endl;
 double delta;
 cin>>delta;
 cout<<"The Result is:"<<endl;
 vector<vector<double> > D;
 vector<vector<double> > L;
 for(int i=0;i<(int)A.size();++i)
 {
 vector<double> z(A.size(),0);
 vector<double> y(A.size(),0);
 for(int j=0;j<(int)A.size();++j)
 {
 if(j<i)
 {
 z[j]=A[j];
 }
 else if(j>i)
 {
 y[j]=A[j];
 }
 }
 D.push_back(y);
 L.push_back(z);
 }
 Print(GS(A,B,D,L,delta));
 system("PAUSE");
}

knight · 发表于 2004-8-4 09:44:10

<

>//松弛迭代（SOR）方法<

>#include <vector>
#include <math.h>
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h><

>using namespace std;ofstream out("out.txt");double operator * (vector<double> a,vector<double> b)
{
 double data=0;
 if(a.size()==b.size())
 {
 vector<double> :: iterator it1=a.begin();
 vector<double> :: iterator it2=b.begin();
 for(;it1!=a.end();++it1,++it2)
 {
 data+=((*it2)*(*it1));
 }
 }
 return data;
}double Delta(vector<double> a,vector<double> b)
{
 double data=0;
 if(a.size()==b.size())
 {
 vector<double> :: iterator it1=a.begin();
 vector<double> :: iterator it2=b.begin();
 for(;(it1!=a.end())&&(it2!=b.end());++it1,++it2)
 {
 if(fabs((*it1)-(*it2))>data)
 {
 data=fabs((*it1)-(*it2));
 }
 }
 }
 return data;
}

vector<double> GS(vector<vector<double> > A,vector<double> B,vector<vector<double> > D,vector<vector<double> > L,double delta,double w)
{
 vector<double> x1(B.size(),0);
 vector<double> x2(B.size(),0);
 do
 {
 x1=x2;
 for(int i=0;i<(int)x2.size();++i)
 {
 double xx=L*x2;
 xx+=D*x1;
 x2=(B-xx)/A;
 }
 for(int i=0;i<(int)x2.size();++i)
 {
 x2=(1-w)*x1+w*x2;
 }
 }while(Delta(x1,x2)>delta);
 return x1;
}void Print(vector<double> a)
{
 vector<double> :: iterator it=a.begin();
 for(;it!=a.end();++it)
 {
 cout<<(*it)<<"\t";
 out<<(*it)<<"\t";
 }
 cout<<endl;
 out<<endl;
}void Print(vector<vector<double> > A)
{
 for(int i=0;i<A.size();++i)
 {
 Print(A);
 }
}

int main()
{
 cout<<"Define the dimension of the matrix:"<<endl;
 int num;
 cin>>num;
 vector<vector<double> > A(num);
 vector<double> B;
 for(int i=0;i<num;++i)
 {
 double bb=0;
 for(int j=0;j<num;++j)
 {
 {
 A.push_back(1/((double)i+(double)j+1));
 bb+=(1/((double)i+(double)j+1));
 }
 }
 B.push_back(bb);
 }
 cout<<"Input the delta:"<<endl;
 double delta;
 cin>>delta;
 cout<<"Input the SOR coefficient:"<<endl;
 double w;
 cin>>w;
 cout<<"The Result is:"<<endl;
 vector<vector<double> > D;
 vector<vector<double> > L;
 for(int i=0;i<(int)A.size();++i)
 {
 vector<double> z(A.size(),0);
 vector<double> y(A.size(),0);
 for(int j=0;j<(int)A.size();++j)
 {
 if(j<i)
 {
 z[j]=A[j];
 }
 else if(j>i)
 {
 y[j]=A[j];
 }
 }
 D.push_back(y);
 L.push_back(z);
 }
 Print(GS(A,B,D,L,delta,w));
 system("PAUSE");
}

knight · 发表于 2004-8-4 09:45:05

以上的C++程序建议用Dev4980编译，其他的编译器本人不敢保证一定编译通过。

knight · 发表于 2004-8-4 09:51:37

<

>不好意思，我忘了，这个程序是我很早以前做得了，大家如果想用还要修改一下，我做得程序是针对Hilbert矩阵来做得（因为当时老师要求用Hilbert矩阵来测试算法），如果要用这个程序算其他的问题请把int main（）里对于矩阵A和向量B的赋值部分修改一下，用你们的A和b来对vector<vector<double>> A和vector<double> B赋值。

freebirdd · 发表于 2004-8-4 17:22:37

[em06]

		自动登录	找回密码
密码			注-册-帐-号