高斯消元法解线性方程组 (高斯-塞得尔算法解方程组)

高斯-塞得尔算法解方程组

{输入方程个数n,再输入n个方程}
{例如方程3x+2y=1,2x+7y=0,则输入:
2
3 2 1
2 7 0
}
{使用了高斯消元法}
program p1052;{解线性方程组}
const maxn=100;
type float=real;
int=integer;
var a:array[1..maxn,1..maxn]of float;
b:array[1..maxn]of float;
d,t,sum:float;
k,l,i,j,m,n,ii,jj,loop:int;
ch:char;
begin
readln(n);
m:=n;
for ii:=1 to n do begin
for jj:=1 to m do read(a[ii][jj]);
read(b[ii]);
end;
k:=1;
d:=0;
l:=0;
while(k<=n)do begin
d:=a[k][k];
l:=k;
for i:=k+1 to n do begin
if(abs(a[i][k])>abs(d)) then begin
d:=a[i][k];
l:=i;
end;
end;
if(l<>k)then begin
for j:=k to n do begin
t:=a[l,j];
a[l,j]:=a[k,j];
a[k,j]:=t;
end;
t:=b[k];
b[k]:=b[l];
b[l]:=t;
end;
for j:=k+1 to n do a[k,j]:=a[k][j]/a[k][k];
b[k]:=b[k]/a[k,k];
for i:=k+1 to n do begin
for j:=k+1 to n do a[i,j]:=a[i,j]-a[i,k]*a[k,j];
j:=1;
b[i]:=b[i]-a[i,k]*b[k];
end;
inc(k);
end;
for i:=n-1 downto 1 do begin
sum:=0;
for j:=i+1 to n do sum:=sum+a[i,j]*b[j];
b[i]:=b[i]-sum;
end;
for loop:=1 to n do begin
if loop=n then ch:=chr(13) else ch:=' ';
write(round(b[loop]),ch);
end;
readln;
readln;
end.

一般认为高斯消元法适用于求解中小规模 (<100元) 的线性方程组.

A 为 m * m 系数矩阵
b 为 m 维常数向量

Code:

#include <cmath>
#include <cstring>
#include <algorithm>
/* Gaussian Elimination for Solving Linear Equations Set */
template <typename T>
inline T * gauss(T * * A, T * b, int m, T eps=1e-15)
{
// [cof con']=[A B']
int p, i, j;
T * * cof=new T *[m];
for (i=0; i<m; i++) {
cof[i]=new T[m];
memcpy(cof[i], A[i], m*sizeof(T));
}
T * con=new T[m];
memcpy(con, b, m*sizeof(T));
T * & x=con;
for (p=0; p<m-1; p++) {
// find major element
int max=p;
for (j=p+1; j<m; j++) {
if (fabs(cof[max][p])<fabs(cof[j][p])) {
max=j;
}
}
if (fabs(cof[max][p])<eps) {
throw "Coefficient matrix is odd";
}
// apply row switching when necessary
if (max!=p) {
std::swap(cof[max], cof[p]);
std::swap(con[max], con[p]);
}
// elimination
for (i=p+1; i<m; i++) {
T rate=cof[i][p]/cof[p][p];
for (j=p; j<m; j++) {
cof[i][j]-=rate*cof[p][j];
}
con[i]-=rate*con[p];
}
}
// back-tracing the unknowns
for (p=m-1; p>0; p--) {
x[p]=con[p]/cof[p][p];
for (j=1; j<=p; j++) {
con[p-j]-=cof[p-j][p]*x[p];
cof[p-j][p]=0.0;
}
}
x[0]=con[0]/cof[0][0];
delete [] cof;
return x;
}



例题, hoj1051 (zju1359) (Beijing 2002 B)

Code:

/**
* Problem: zju1359
* Source: ACM/ICPC Asia Regional (Beijing), 2002, B. Random Walk
* Author: LIU Yu, <pineapple.liu@gmail.com>
* Status: AC
*/

#include <cstdio>
#include <cmath>
#include <cctype>
#include <cstring>
#include <string>
#include <map>
#include <vector>
#include <algorithm>

template <typename T>
inline T * gauss(T * * A, T * b, int m, T eps=1e-15)
{
// [cof con']=[A B']
int p, i, j;
T * * cof=new T *[m];
for (i=0; i<m; i++) {
cof[i]=new T[m];
memcpy(cof[i], A[i], m*sizeof(T));
}
T * con=new T[m];
memcpy(con, b, m*sizeof(T));
T * & x=con;
for (p=0; p<m-1; p++) {
// find major element
int max=p;
for (j=p+1; j<m; j++) {
if (fabs(cof[max][p])<fabs(cof[j][p])) {
max=j;
}
}
if (fabs(cof[max][p])<eps) {
throw "Coefficient matrix is odd";
}
// apply row switching when necessary
if (max!=p) {
std::swap(cof[max], cof[p]);
std::swap(con[max], con[p]);
}
// elimination
for (i=p+1; i<m; i++) {
T rate=cof[i][p]/cof[p][p];
for (j=p; j<m; j++) {
cof[i][j]-=rate*cof[p][j];
}
con[i]-=rate*con[p];
}
}
// back-tracing the unknowns
for (p=m-1; p>0; p--) {
x[p]=con[p]/cof[p][p];
for (j=1; j<=p; j++) {
con[p-j]-=cof[p-j][p]*x[p];
cof[p-j][p]=0.0;
}
}
x[0]=con[0]/cof[0][0];
delete [] cof;
return x;
}

typedef struct
{
enum
{
IF_FAR, IF_NEAR, NOP,
} type;
union
{
struct
{
double prob;
int target;
} if_data;
};
} command_t;

void emulate(int index, command_t * proc, int loc,
int eip, double factor, double * a, double& b)
{
if (eip >= loc) return;
if (factor < 1e-10) return;

// time consumed by current instruction
b += factor;

command_t & cmd = proc[eip];
switch (cmd.type)
{
case command_t::IF_FAR:
a[cmd.if_data.target] -= factor * cmd.if_data.prob;
break;
case command_t::IF_NEAR:
if (cmd.if_data.target == 1)
{
a[index] -= factor * cmd.if_data.prob;
}
else
{
emulate(index, proc, loc, cmd.if_data.target - 1,
factor * cmd.if_data.prob, a, b);
}
factor *= 1.0 - cmd.if_data.prob;
break;
case command_t::NOP:
break;
default:
return;
}
emulate(index, proc, loc, eip + 1, factor, a, b);
}

#define PROG_START "PROG_START"
#define PROG_END "PROG_END"
#define PROC_START "PROC"
#define PROC_END "END"
#define CMD_IF "IF"
#define CMD_NOP "NOP"
#define REQUEST_END "REQUEST_END"

int main(void)
{
using namespace std;

map<string, int> dict;

int N = 0;
command_t proc[100][100];
int loc[100] = {0};

char line[1024];
while (true)
{
gets(line);
if (strncmp(line, PROG_START, strlen(PROG_START)) == 0)
{
continue;
}
if (strncmp(line, PROG_END, strlen(PROG_END)) == 0)
{
break;
}
if (strncmp(line, PROC_START, strlen(PROC_START)) == 0)
{
char name[1024];
sscanf(line, PROC_START " %s", name);
dict[name] = N;
while (true)
{
gets(line);
if (strncmp(line, PROC_END, strlen(PROC_END)) == 0)
{
break;
}
command_t cmd;
if (strncmp(line, CMD_NOP, strlen(CMD_NOP)) == 0)
{
cmd.type = command_t::NOP;
}
else if (strncmp(line, CMD_IF, strlen(CMD_IF)) == 0)
{
char cmp;
double prob;
char action[1024];
char target[1024];
sscanf(line, CMD_IF " x%c%lf %s %[A-Za-z0-9]s",
&cmp, &prob, action, target);
cmd.type = (action[0] == 'P')?(command_t::IF_FAR):
(command_t::IF_NEAR);
cmd.if_data.prob = (cmp == '>')?(1.0 - prob):(prob);
cmd.if_data.target = (isalpha(target[0]))?(dict[target]):
(atoi(target));
}
proc[N][loc[N]++] = cmd;
}
N++;
}
}

double * * A = new double * [N];
double * B = new double[N];
for (int i = 0; i < N; i++)
{
A[i] = new double[N];
memset(A[i], 0, N * sizeof(double));
A[i][i] += 1.0;
B[i] = 0.0;
emulate(i, proc[i], loc[i], 0, 1.0, A[i], B[i]);
}

double * X = gauss(A, B, N, 1e-4);

while (true)
{
gets(line);
if (strncmp(line, REQUEST_END, strlen(REQUEST_END)) == 0)
{
break;
}
char name[1024];
sscanf(line, "%[A-Za-z0-9]s", name);
printf("%.3lf\n", X[dict[name]]);
}

delete [] X;
delete [] B;
for (int j = 0; j < N; j++)
{
delete [] A[j];
}
delete [] A;

return 0;
}





第22行和第26行cof[p][max]、cof[p][j]应该换成cof[max][p], cof[j][p]。
容易找到 让原模版出错的数据
0x+1y+0z = 1
0x+0y+1z = 1
1x+0y+0z = 1