雅克比迭代法……

1、实验课程名称：数值分析 实验项目名称 运用插值法 实验成绩 实验者 江骏 专业班级 软件0803 组别 同组者 实验日期 年 月 日 第一部分：实验分析与设计 一、实验内容描述 研究用Jacobi迭代法与 Gauss-Seidel迭代法解下列方程组Ax=b的收敛性，通过上机计算，验证分析是 否正确，并观察右端项对迭代收敛是否有影响，比较两法的收敛速度 松弛因子对超松弛因子迭代法收敛速度的影响，要求对不同的阶数进行迭代 观察欧拉显式方法的收敛性 观察欧拉隐式方法的收敛性 写出实验报告 二、实验基本原理与设计 Jacobi迭代法：刀1=1玄耳xj=b(i=1,2,n) Gauss-Seidel 迭

2、代法：x i(K+1) =b i-刀 i-1 j=1 ajXj(k+1)-刀爲+1为(k) /a j 超松弛因子迭代法：Dx (k+1) =Dx (k) +Ex (k+1) +Fx (k) +b-Dx (k) X(k+1) =x(k)+D-1R 三角分解法：利用三角矩阵将方程组化解为两个方程组，从而简化运算。 欧拉显式方法：在任意节点tn+1=to+(n+1)h处，u(t n+1)的近似值由Euler公式给出： Un+1=Un+hf(t n,u n) 欧拉隐式方法：改进的 Euler公式y n+1=yn+hf(x n,y n) h yn+1=yn+ f(x n,y n)+f(x n+1, y

3、n+1) 2 三、主要仪器设备及耗材 Win dowsXp VC+6.0 第二部分：实验调试与结果分析 一、调试过程 Jacobi迭代法 #in clude #defi ne N 10 float ABS(float,float); int main (void) int i, j, n; float aNN, bN; float xN, yN; float e,total; printf (Please in put the dime nsion:”)； scanf (%d, printf (Please in put the Coefficie nt Matrix:); for (i =

4、0; i n; i+) for (j = 0; j n; j+) scanf (%f, printf (Please in put the Vector:); for (i = 0; i n; i+) scanf (”f, printf (Please in put the in itial vector:); for (i = 0; i n; i+) scanf (%f, do for (i = 0; i n; i+) total = 0.0f; for (j = 0; j n; j+) if (i != j) total += aij * xj; yi = (bi - total) / a

5、ii; e=0.0f; for (j = 0; j n; j+) e = e + ABS(xj, yj); prin tf(%fn,e); for (i = 0; i 0.0001); for (i = 0; i n; i+) printf (x%d=%ft, i, xi); printf (n); return 0; float ABS(float x,float y) int total; if(xy) total=y-x; else total=x-y; return total; 结果 cQ Docuaents and SettincsXAd M=ur)上 u_d 宀 二一X-J_%=

6、u_d 宀 sx * _m=p = T+ov 一O.II)O4 )(puno)-ZM 宀 宀 兰曰殳日三亨-0 JlrH _m二p (+OVOUDO4 兰三e二二A 曰6 T+ov 一O.II)O4 :=?lnH puno二 u- -56-55P Qqnop 二寸 Ze-) HSA Qqnop 二寸-Lew-寸-L-VZ9 H55B Qqnop 超松弛因子迭代法 #in elude #in elude using n amespaee std; float *on e_array_malloc(i nt n); float *two_array_malloc(i nt m,i nt n); f

7、loat matrix_category(float* x,i nt n); int mai n() con st i nt MAX=100; int n,i,j,k; float* a; float* x_0; float* x_k; float precisi on; float w; cout precisi on; coutendl输入系数矩阵的阶数，N： cinn; a=two_array_malloc (n,n+1); coutendl输入增广矩阵的各值：n; for(i=0;i n;i+) for(j=0;j aij; x 0=on e array malloc (n); cou

8、tendl 输入初始向量：n; for(i=0;ix_Oi; x_k=on e_array_malloc( n); cout输入松弛因子 w (1w w; float temp; for(k=0;kMAX;k+) for(i=0;i n;i+) temp=0; for(j=0;ji;j+) temp=temp+aij*x_kj; x_ki=ai n-temp; temp=0; for(j=i+1;j n;j+) temp=temp+aij*x_0j; x_ki=(x_ki-temp)/aii; x_ki=(1-w)*x_0i+w*x_ki; for(i=0;i n;i+) x_0i=x_ki-

9、x_0i; if(matrix_category(x_ 0,n) vprecisi on) break; else for(i=0;i n;i+) x_0i=x_ki; if(MAX=k) cout迭代不收敛n; cout迭代次数为：ke ndl; cout解向量为：n; for(i=0;i n;i+) coutxi: x_kie ndl; return 0; float *on e_array_malloc(i nt n) float *a; a=(float *)malloc(sizeof(float)* n); return a; float *two_array_malloc(i nt

10、 m,i nt n) float *a; int i; a=(float *)malloc(m*sizeof(float *); for (i=0;im;i+) ai=(float *)malloc( n*sizeof(float); return a; float matrix_category(float* x,i nt n) int i; float temp=0; for(i=0;i n;i+) temp=temp+fabs(xi); return temp; nllx 结果： 流程图: 三角分解法 #in clude using n amespace std; int mai n()

11、 const int N=100; static double aNN,bN; int i,j,k,num,p; double m,t,q; cout请输入矩阵阶数： cinnum; / 输入 aij,bi cout now in put the matrix aij,i,j=1 num:e ndl; for(i=1;i=nu m;i+) for(j=1;j=nu m;j+) coutaijaij; cout now in put the matrix bi,i=1 num :e ndl; for(i=1;i=nu m;i+) coutbi bi; t=0; for(i=1;i=nu m;i+

12、) m=0; for(j=i;j=nu m+1;j+) for(k=1;k=i-1;k+) m=m+aik*akj; aij= aij-m; for(j=i+1;j=nu m;j+) for(k=1;k=1;k-) for(j=k+1;j=nu m;j+) q=q+akj*aj num+1; ak num+1=(ak nu m+1-q)/akk; cout方程组的解为：; for(i=1;i=nu m;i+) coutxi=ai num+1ents and SettingsAdinistratorDebugsanjiao.exe* 歸1翩数:M nou input the matrix aij

13、,i,j=l.3 : klll=6 2 1 1 4 -2 -3 1 4 kl2-al3-a2l-a22-a23-a3l-a32-a33-now input ths ma tpix bi,i=1.3 : bll=-3 2 4 k2-b3-方程组的解为：xlJ=0.0697509 L21=0.0520741 Lc3 3=0.165049 Press any key to continue 流程图: 开始 矩阵阶数 系数矩阵 for(k=1;k=i-1;k+) m=m+aik*akj; aij= aij-m; ak nu m+1=(ak num+1-q)/akk 结束 欧拉显式方法 #in elud

14、e #in elude #in elude double f(double x,double y) return x*pow(y,1.0/3); int mai n() int i; double x,y,y0=1,dx=0.1; double xx11; double euler11,euler_211; double temp; double f(double x,double y); for (i=0;i11;i+) xxi=i*dx; euler0=y0; for (i=1,x=0;i11;i+,x+=dx) euleri=euleri-1+dx*f(x,euleri-1); eule

15、r_20=y0; for (i=1,x=0;i11;i+,x+=dx) temp=euler_2i-1+dx*f(x,euler_2i-1); euler_2i=euler_2i-1+dx*(f(x,euler_2i-1)+f(x+dx,temp)/2; for (i=0,x=0;i11;i+,x+=dx) prin tf(x=%lfteluer=%lfteuler_2=%lftaccu=%lfn,x,euleri,euler_2i,pow(1+x*x,1.0/3); getch(); 流程图： 欧拉隐式方法 #in clude #in clude #defi ne f(x,y)(x+y) i

16、nt mai n() int m; int i; double a,b,y0; double xn,yn,xn1,yn1,yn 1b; double h; prin tf(nln put the begi n and end of x:); scan f(%lf%lf, prin tf(I nput the y value at %f:,a); scan f(%lf, prin tf(I nput m valuedivide(%f,%f):,a,b); scan f(%d, if(m=0) pri ntf(Please in put a nu mber larger tha n 1.n); r

17、eturn(1); h=(b-a)/m; xn=a;yn=yO; for(i=1;i=m;i+) xn 1=x n+h; yn 1b=yn+h*f(x n,yn); yn 1=yn+h/2*(f(x n,y n)+f(x n1,y n1b); prin tf(x%d=%f,y%d=%fn,i,xxn 1,i,y n1); xn=xn1; yn=yn1; return(0); 流程图： 二、实验结果及分析 雅克比迭代法的全部分量都是间接利用，由于新值比旧值好，Gauss-Seidel迭代法直接对算出 的分量加以使用，这样使得迭代的收敛情况有所改善，而超松弛法是在Gauss-Seidel迭代法的 基础上，应