下面是小编为大家整理的雅克比迭代法(2022年),供大家参考。
实验项目名称
运用插值法 实验成绩
实验者
江骏 专业班级
软件 0803 组别
同组者
第一部分:实验分析与设计
一、实验内容描述 实验日期
年 月 日
⑴ 研究用
Jacobi
迭代法与
Gauss- - Seidel
迭代法解下列方程组
Ax=b
的收敛性,通过上机计算,验证分析是否正确,并观察右端项对迭代收敛是否有影响,比较两法的收敛速度
⑵ 松弛因子对超松弛因子迭代法收敛速度的影响,要求对不同的阶数进行迭代
⑶ 观察欧拉显式方法的收敛性
⑷ 观察欧拉隐式方法的收敛性
⑸ 写出实验报告
二、实验基本原理与设计
①Jacobi
迭代法:
∑ n n
a a
x x
=b
(i=1,2,…,n)
j=1
ij
j j
i i
②Gauss- - Seidel
迭代法 :x x
(K+1)
=[b
- - ∑ i i- -1 1
a x
(k+1)
- - ∑ n n
a x (k)
]/a
i i
i i
j=1
ij
j j
j=i+1
ij
j j
ij
③ 超松弛因子迭代法:
Dx (k+1)
=Dx (k)
+Ex
(k+1)
+Fx (k)
+b- - Dx
(k)
X X (k+1) =x (k) +D - -1 1
R R
④ 三角分解法:利用三角矩阵将方程组化解为两个方程组,从而简化运算。
⑤ 欧拉显式方法:在任意节 点 t t
=t
+(n+1)h
处, u(t
) ) 的近似值由
Euler
公式给出:
n+1
0 0
n+1
u u
=u
+hf(t
,u
) )
n+1
n n
n n
n n
⑥ 欧拉隐式方法:改进 的 Euler
公式
y =y
+hf(x
,y
) )
n+1
n n
n n
n n
h y y
=y
+ +
[f(x
,y
)+f(x
, ,
y )]
n+1
n n
2 n n
n n
n+1
n+1
三、主要仪器设备及耗材
WindowsXp VC++
第二部分:实验调试与结果分析
一、调试过程
①Jacobi 迭代法#include <> #define N 10
float ABS(float,float); int
main (void)
{ {
int i, j, n;
float a[N][N], b[N];
float x[N], y[N]; float e,total;
printf ("Please input the dimension:"); scanf ("%d", &n);
printf ("Please input the Coefficient Matrix: "); for (i = 0; i < n; i++)
for (j = 0; j < n; j++) scanf ("%f", &a[i][j]);
printf ("Please input the Vector: "); for (i = 0; i < n; i++)
scanf ("%f", &b[i]);
printf ("Please input the initial vector: "); for (i = 0; i < n; i++)
scanf ("%f", &x[i]);
do
{ {
for (i = 0; i < n; i++)
{ {
total = ;
for (j = 0; j < n; j++)
{ {
if (i != j)
total += a[i][j] * x[j];
} }
y[i] = (b[i] - -
total) / a[i][i];
} }
e=;
for (j = 0; j < n; j++)
e = e + ABS(x[j], y[j]);
printf("%f\ \ n",e);
for (i = 0; i < n; i++)
x[i] = y[i];
} }
while (e > ;
for (i = 0; i < n; i++)
printf ("x%d=%f\ \ t", i, x[i]); printf ("\ \ n");
return 0;
} }
float ABS(float x,float y)
{ {
int total;
if(x<y) total=y- - x; else total=x- - y; return total;
} }
结 果 :
流程图:
开始
阶数 n,系数矩阵,初始向量
for(i=0;i<n;i+
Y
i=j
N
total
+=
e e
= =
e e
+ +
结束
②Gauss- - Seidel
迭代法#include<>
int
main()
{ {
double x[3] = {0, 0, 0};
double a[3][3] = {6,2,- - 1,1,4,- - 2,- - 3,1,4};
double y[3] = {- - 3,2,4};
double d[3][3],g[3]; int round = 5, i,j; for (i=0; i<3; ++i)
{ {
g[i] = y[i] / a[i][i]; for (j=0; j<3; ++j)
{ {
d[i][j] = i==j ? 0 : - - a[i][j]/a[i][i];
} }
} }
while (round -- ) { for (i=0; i<3;
++i)
{ {
x[i] = g[i];
for (j=0; j<3; ++j)
{ {
x[i] += d[i][j] * x[j];
} }
printf("%lf", x[i]);
} }
printf("\ \ n");
} }
} }
结 果 :
流程图:
开始
系数矩阵
for
(i=0;
i<3;
++i) d[i][j]
= =
i==j
? ?
0 0
: :
for
(j=0;
j<3;
++j)
结束
③ 超松弛因子迭代法#include <iostream> #include <cmath> using namespace std;
float *one_array_malloc(int n);
float **two_array_malloc(int m,int n); float matrix_category(float* x,int n); int main(){ const int MAX=100;
int n,i,j,k; float** a; float* x_0; float* x_k;
float
precision;
float w; cout<<" 输
入
精
度
e:"; cin>>precision;
cout<<endl<<" 输入系数矩阵的 阶数,N N :
"; cin>>n; a=two_array_malloc(n,n+1); cout<<endl<<" 输入增广矩阵的各值:\ \ n"; for(i=0;i<n;i++)
{ {
for(j=0;j<n+1;j++)
{ {
cin>>a[i][j];
} }
} }
x_0=one_array_malloc(n); cout<<endl<<" 输入初始向量: :\ \ n"; for(i=0;i<n;i++)
{ {
cin>>x_0[i]; } x_k=one_array_malloc(n);
cout<<" 输入松弛因子w w
(1<w<2) :\ \ n"; cin>>w;
float temp; for(k=0;k<MAX;k++)
{ {
for(i=0;i<n;i++)
{ {
temp=0; for(j=0;j<i;j++)
{ {
temp=temp+a[i][j]*x_k[j];
} }
x_k[i]=a[i][n]- - temp; temp=0; for(j=i+1;j<n;j++)
{ {
temp=temp+a[i][j]*x_0[j];
} }
x_k[i]=(x_k[i]- - temp)/a[i][i];
x_k[i]=(1- - w)*x_0[i]+w*x_k[i];
} }
for(i=0;i<n;i++)
{ {
x_0[i]=x_k[i]- - x_0[i];
} if(matrix_category(x_0,n)<precision)
{ {
break;
} }
else
{ {
for(i=0;i<n;i++)
{ {
x_0[i]=x_k[i];
} }
} }
} }
if(MAX==k)
{ {
cout<<" 迭代不收敛\ \ n";
} cout<<" 迭代次数为:
"<<k<<endl; cout<<" 解
向
量
为
:\ \ n"; for(i=0;i<n;i++)
{ {
cout<<"x"<<i<<":
"<<x_k[i]<<endl;
} return
0;
} }
float *one_array_malloc(int n)
{ {
float
*a;
a=(float *)malloc(sizeof(float)*n); return a;
} }
float **two_array_malloc(int m,int n)
{ {
float **a; int
i;
a=(float **)malloc(m*sizeof(float *)); for (i=0;i<m;i++)
{ {
a[i]=(float
*)malloc(n*sizeof(float));
} }
return a;
} }
float matrix_category(float* x,int
n)
{ int i; float temp=0;
for(i=0;i<n;i++)
{ {
temp=temp+fabs(x[i]);
} }
return temp;
} }
结果:
流程图:
开始
系数矩阵阶, 精度,增广阵,初始
for(k=0;k<MAX;
Y
matrix_category(x
N
for(i=0;i<n;i++)
a=(float
*)malloc(sizeof
结束
④ 三角分解法#include<iostream> using namespace std; int main()
{ {
const int N=100;
static
double
a[N][N],b[N];
int i,j,k,num,p; double m,t,q;
cout<<" 请输入矩阵阶数:
"; cin>>num;
."<<num<<" :"<<endl; for(i=1;i<=num;i++)
{ {
for(j=1;j<=num;j++)
{ {
cout<<"a["<<i<<"]["<<j<<"]=";
cin>>a[i][j];
} }
} }
cout<<"now input the matrix b[i],i=1..."<<num<<" :"<<endl; for(i=1;i<=num;i++)
{ {
cout<<"b["<<i<<"]="; cin>>b[i];
} }
t=0;
for(i=1;i<=num;i++)
{ {
m=0;
for(j=i;j<=num+1;j++)
{ {
for(k=1;k<=i- - 1;k++)
m=m+a[i][k]*a[k][j];
a[i][j]=
a[i][j]- - m;
} }
for(j=i+1;j<=num;j++)
{ {
for(k=1;k<=i- - 1;k++) t=t+a[j][k]*a[k][i];
a[j][i]=(a[j][i]- - t)/a[i][i];
} }
} }
a[num][num+1]=a[num][num+1]/a[num][num]; q=0;
for(k=num- - 1;k>=1;k --) )
{ {
for(j=k+1;j<=num;j++)
q=q+a[k][j]*a[j][num+1];
a[k][num+1]=(a[k][num+1]- - q)/a[k][k];
} }
cout<<" 方
程
组
的
解
为
:
"; for(i=1;i<=num;i++) cout<<"x["<<i<<"]="<<a[i][num+1]<<endl; return 0;
} }
结 果 :
流程图:
开始
矩阵阶数
for(k=1;k<=i- - 1;k++) m=m+a[i][k]*a[k][j];
a[k][num+1]=(a[k][
结束
⑤ 欧拉显式方法#include <> #include <> #include <>
double f(double x,double y)
{ {
return x*pow(y,3);
} }
int main()
{ {
int i;
double x,y,y0=1,dx=; double xx[11];
double euler[11],euler_2[11];
double temp;
double f(double x,double y); for (i=0;i<11;i++) xx[i]=i*dx;
euler[0]=y0;
for (i=1,x=0;i<11;i++,x+=dx) euler[i]=euler[i- - 1]+dx*f(x,euler[i- - 1]); euler_2[0]=y0;
for (i=1,x=0;i<11;i++,x+=dx)
{ {
temp=euler_2[i- - 1]+dx*f(x,euler_2[i- - 1]);
euler_2[i]=euler_2[i- - 1]+dx*(f(x,euler_2[i- - 1])+f(x+dx,temp))/2;
} }
for (i=0,x=0;i<11;i++,x+=dx) printf("x=%lf\ \ teluer=%lf\ \ teuler_2=%lf\ \ taccu=%lf\ \ n",x,euler[i],euler_2[i],pow(1+x*x,3)); getch();
} }
结 果 :
流程图:
⑥ 欧拉隐式方法
#include<> #include<>
#define f(x,y)(x+y ) int
main()
{int m; int i;
double a,b,y0;
double xn,yn,xn1,yn1,yn1b; double h;
printf("\ \ nInput the begin and end of x:"); scanf("%lf%lf",&a,&b);
printf("Input the y value at %f:",a); scanf("%lf",&y0);
printf("Input m value[divide(%f,%f)]:",a,b); scanf("%d",&m);
if(m<=0)
{printf("Please input a number larger than 1.\ \ n"); return(1);
} }
h=(b- - a)/m; xn=a;yn=y0; for(i=1;i<=m;i++)
{xn1=xn+h; yn1b=yn+h*f(xn,yn);
yn1=yn+h/2*(f(xn,yn)+f(xn1,yn1b)); printf("x%d=% f,y%d=%f\ \ n",i,xn1,i,yn1); xn=xn1;
yn=yn1;
}return(0);
} }
结 果 :
流程图:
二、实验结果及分析 雅克比迭代法的全部分量都是间接利用,由于新值比旧值好,Gauss-Seidel 迭代法直接对算出
的分量加以使用,这样使得迭代的收敛情况有所改善,而超松弛法是在 Gauss-Seidel 迭代法的 基础上,应用加速收敛思想而构造的一种方法;欧拉公式和梯形公式在计算上有明显的区别, 欧拉公式的特点是可以有 U 直接计算 U ,也就是显式的,而梯形公式右端也有 U ,必须通过解 N n+1 n+1
方程才能求得,这类是隐式的。从算法的角度看显式远比隐式方便,隐式通常用迭代法求解。
三、实验小节、建议及体会 任何一种方法都有它的优点和局限性,只有找到一种方法能够解决我们的问题,又尽量减少出错,减少 计算量就是符合这个问题的好方法。
良好的编程习惯有助于检查错误,也可以增强程序的可阅读性。