- 技术人生
- 高斯消元法
高斯消元法
方法1：
#include
"stdio.h"
#include
"stdlib.h"
void
RKT
(
t
,
y
,
n
,
h
,
k
,
z
)
int
n
;
/*
微分方程组中方程的个数，也是未知函数的个数
*/
int
k
;
/*
积分的步数
(
包括起始点这一步
)*/
double
t
;
/*
积分的起始点
t0*/
double
h
;
/*
积分的步长
*/
double
y
[];
/*
存放
n
个未知函数在起始点
t
处的函数值
,
返回时
,
其初值在二维数组
z
的第零列中
*/
double
z
[];
/*
二维数组
,
体积为
n x k.
返回
k
个积分点上的
n
个未知函数值
*/
{
extern
void
Func
();
/*
声明要求解的微分方程组
*/
int
i
,
j
,
l
;
double
a
[4],*
b
,*
d
;
b
=
malloc
(
n
*
sizeof
(
double
));
/*
分配存储空间
*/
if
(
b
==
NULL
)
{
printf
(
"
内存分配失败
/n"
);
exit
(1);
}
d
=
malloc
(
n
*
sizeof
(
double
));
/*
分配存储空间
*/
if
(
d
==
NULL
)
{
printf
(
"
内存分配失败
/n"
);
exit
(1);
}
/*
后面应用
RK4
公式中用到的系数
*/
a
[0]=
h
/2.0;
a
[1]=
h
/2.0;
a
[2]=
h
;
a
[3]=
h
;
for
(
i
=0;
i
<=
n
-1;
i
++)
z
[
i
*
k
]=
y
[
i
];
/*
将初值赋给数组
z
的相应位置
*/
for
(
l
=1;
l
<=
k
-1;
l
++)
{
Func
(
y
,
d
);
for
(
i
=0;
i
<=
n
-1;
i
++)
b
[
i
]=
y
[
i
];
for
(
j
=0;
j
<=2;
j
++)
{
for
(
i
=0;
i
<=
n
-1;
i
++)
{
y
[
i
]=
z
[
i
*
k
+
l
-1]+
a
[
j
]*
d
[
i
];
b
[
i
]=
b
[
i
]+
a
[
j
+1]*
d
[
i
]/3.0;
}
Func
(
y
,
d
);
}
for
(
i
=0;
i
<=
n
-1;
i
++)
y
[
i
]=
b
[
i
]+
h
*
d
[
i
]/6.0;
for
(
i
=0;
i
<=
n
-1;
i
++)
z
[
i
*
k
+
l
]=
y
[
i
];
t
=
t
+
h
;
}
free
(
b
);
/*
释放存储空间
*/
free
(
d
);
/*
释放存储空间
*/
return
;
}
main
()
{
int
i
,
j
;
double
t
,
h
,
y
[3],
z
[3][11];
y
[0]=-1.0;
y
[1]=0.0;
y
[2]=1.0;
t
=0.0;
h
=0.01;
RKT
(
t
,
y
,3,
h
,11,
z
);
printf
(
"/n"
);
for
(
i
=0;
i
<=10;
i
++)
/*
打印输出结果
*/
{
t
=
i
*
h
;
printf
(
"t=%5.2f/t "
,
t
);
for
(
j
=0;
j
<=2;
j
++)
printf
(
"y(%d)=%e "
,
j
,
z
[
j
][
i
]);
printf
(
"/n"
);
}
}
void
Func
(
y
,
d
)
double
y
[],
d
[];
{
d
[0]=
y
[1];
/*y0'=y1*/
d
[1]=-
y
[0];
/*y1'=y0*/
d
[2]=-
y
[2];
/*y2'=y2*/
return
;
}
高斯消元法：
#include
<
stdio
.h>
#include
<stdlib.h>
#include
<
malloc
.h>
#include
<math.h>
int
GS
(
int
,
double
**,
double
*,
double
);
double
**
TwoArrayAlloc
(
int
,
int
);
void
TwoArrayFree
(
double
**);
void
main
()
{
int
i
,
n
;
double
ep
,**
a
,*
b
;
n
= 3;
ep
= 1e-4;
a
=
TwoArrayAlloc
(
n
,
n
);
b
= (
double
*)
calloc
(
n
,
sizeof
(
double
));
if
(
b
==
NULL
)
{
printf
(
"
内存分配失败
/n"
);
exit
(1);
}
a
[0][0]= 2;
a
[0][1]= 6;
a
[0][2]=-1;
a
[1][0]= 5;
a
[1][1]=-1;
a
[1][2]= 2;
a
[2][0]=-3;
a
[2][1]=-4;
a
[2][2]= 1;
b
[0] = -12;
b
[1] = 29;
b
[2] = 5;
if
(!
GS
(
n
,
a
,
b
,
ep
))
{
printf
(
"
不可以用高斯消去法求解
/n"
);
exit
(0);
}
printf
(
"
该方程组的解为：
/n"
);
for
(
i
=0;
i
<3;
i
++)
printf
(
"x%d = %.2f/n"
,
i
,
b
[
i
]);
TwoArrayFree
(
a
);
free
(
b
);
}
int
GS
(
n
,
a
,
b
,
ep
)
int
n
;
double
**
a
;
double
*
b
;
double
ep
;
{
int
i
,
j
,
k
,
l
;
double
t
;
for
(
k
=1;
k
<=
n
;
k
++)
{
for
(
l
=
k
;
l
<=
n
;
l
++)
if
(
fabs
(
a
[
l
-1][
k
-1])>
ep
)
break
;
else
if
(
l
==
n
)
return
(0);
if
(
l
!=
k
)
{
for
(
j
=
k
;
j
<=
n
;
j
++)
{
t
=
a
[
k
-1][
j
-1];
a
[
k
-1][
j
-1]=
a
[
l
-1][
j
-1];
a
[
l
-1][
j
-1]=
t
;
}
t
=
b
[
k
-1];
b
[
k
-1]=
b
[
l
-1];
b
[
l
-1]=
t
;
}
t
=1/
a
[
k
-1][
k
-1];
for
(
j
=
k
+1;
j
<=
n
;
j
++)
a
[
k
-1][
j
-1]=
t
*
a
[
k
-1][
j
-1];
b
[
k
-1]*=
t
;
for
(
i
=
k
+1;
i
<=
n
;
i
++)
{
for
(
j
=
k
+1;
j
<=
n
;
j
++)
a
[
i
-1][
j
-1]-=
a
[
i
-1][
k
-1]*
a
[
k
-1][
j
-1];
b
[
i
-1]-=
a
[
i
-1][
k
-1]*
b
[
k
-1];
}
}
for
(
i
=
n
-1;
i
>=1;
i
--)
for
(
j
=
i
+1;
j
<=
n
;
j
++)
b
[
i
-1]-=
a
[
i
-1][
j
-1]*
b
[
j
-1];
return
(1);
}
double
**
TwoArrayAlloc
(
int
r
,
int
c
)
{
double
*
x
,**
y
;
int
n
;
x
=(
double
*)
calloc
(
r
*
c
,
sizeof
(
double
));
y
=(
double
**)
calloc
(
r
,
sizeof
(
double
*));
for
(
n
=0;
n
<=
r
-1;++
n
)
y
[
n
]=&
x
[
c
*
n
];
return
(
y
);
}
void
TwoArrayFree
(
double
**
x
)
{
free
(
x
[0]);
free
(
x
);
}
