計算物理-雅克比迭代法解方程組

你想擺脫三元四元甚至五元一次方程組的計算嗎？你想讓自己無論遇到什麼樣的一次方程組都能按照通用的套路把精確的解求出來嗎？來吧，看我用程序寫一個通用的解決方程組求解的問題，無論你是多少組方程，我都能求出來！求出來的速度僅僅取決於你的電腦處理器的速度。

迭代法是一種不斷套用一個迭代公式 , 逐步逼近方程組的解的方法 , 特別適合於大型稀疏矩陣 . 常用的迭代法有雅可比 (Jacobi) 迭代法 ( 又稱簡單迭代法 ) 和賽德爾 (Seidel) 迭代法 .

設線性代數方程組的一般形式為

逐一變數分離得到：

則最後代出來的表達式為：

原來的方程組是數組{A1.....An}=B1，其中An和Bn都已知。我們需要求解的是數組{X1,X2,X3....Xn}，由新的表達式可以得出。我們需要求解的數組由任意一個初始數組設定值按照要求的方程組帶入之後再將計算得到的新的值重新帶入原方程組求解，重複以上過程之後，第i次計算得到的結果Xi和第i+1次計算得到的結果相減，即X(i)-X(i+1)<a，其中a為我們設定的精度值，這樣可以得到符合精度要求的方程組的數值解。

但是在求解數值解之前需要先確定迭代的收斂和發散情況，如果如下極限存在，則迭代是收斂的，也就是說可無數次迭代之後數值趨於某個穩定值不變。如果不存在，這說明迭代過程是發散的。

代碼段：

#include "StdAfx.h"//雅克比迭代法輪子1n#include "stdlib.h"n#include "stdio.h"n#include "conio.h"n#include "string.h"n#include "math.h"n#define N 100nfloat Table(int n,float a[N][N],float b[N]) n{ntint i,j;ntfloat c[N][N];ntprintf("請輸入矩陣的行!n");n for(i=0;i<n;i++)nt{nttprintf("請輸入行 %d:",i+1);nttfor(j=0;j<n;j++)ntttscanf("%f",&a[i][j]); nt}ntprintf("請輸入向量b:");ntfor(i=0;i<n;i++)nttscanf("%f",&b[i]); n for(i=0;i<n;i++)n for(j=0;j<n;j++)n {n if(i==j)n {n c[i][j]=0;n continue;n }n c[i][j]=-a[i][j]/a[i][i]; n }n printf("n矩陣A和向量b:n");ntfor(i=0;i<n;i++)nt{nttfor(j=0;j<n;j++)ntttprintf("%10.5f",a[i][j]);n printf("%10.5f",b[i]);ntt printf("n");nt}nntprintf("n高斯-賽德爾迭代的方案:n");n for(i=0;i<n;i++)n {n for(j=0;j<n;j++)n printf("%10.5f",c[i][j]);n printf("%10.5f",b[i]/a[i][i]);n printf("n");n }n}nfloat init_vec(int n,float x[N])n{n int i;n printf("請輸入初始迭代向量x:");n for(i=0;i<n;i++)n scanf("%f",&x[i]); n printf("n初始迭代向量x:n");n for(i=0;i<n;i++)n printf("%10.5f",x[i]);n printf("n");n}nfloat jacobi(int n,float a[N][N],float b[N],float x[N]) n{n int i,j,k;n float tmp,x2[N];n for(k=0;;k++)n {n for(i=0;i<n;i++)n x2[i]=x[i];n for(i=0;i<n;i++)nttt{n tmp=0.0;n for(j=0;j<n;j++)n {n if(j==i) continue;n tmp+=a[i][j]*x2[j];n }n x[i]=(b[i]-tmp)/a[i][i]; n }n for(i=0,j=0;i<n;i++)n if(fabs(x2[i]-x[i])<0.00001) j++; n if(j==n) n {n printf("n這個高斯-賽德爾迭代方案收斂!n");n printf("迭代次數:%d",k+1); n break;n }n if(k==499) n {n printf("n這雅可比迭代計劃可能不收斂!");n break; n }n }n printf("n結果為：n");n for(i=0;i<n;i++)n printf("%12.7f",x[i]); n}nint main()n{n int n;n float x[N],a[N][N],b[N];n printf("請輸入n:");n scanf("%d",&n); n Table(n,a,b); n init_vec(n,x); n jacobi(n,a,b,x); n getch();n}n

#include "stdafx.h"//雅克比迭代輪子2n#include<stdio.h>n#include<math.h>nint main(void)n{n double A[5][5] = {{28,-3,0,0,0},n {-3,38,-10,0,-5},n {0,-10,25,-15,0},n {0,0,-15,45,0},n {0,-5,0,0,30}};n double b[5] = {10,0,0,0,0};n double x[5] = {0}; //第k+1次迭代的結果n double xx[5] = {0}; //第k次迭代的結果n int size = 5;n int Max = 100; //最大迭代次數n double residual = 0.0; //n double sum = 0.0; n double dis = 0.0;n double dif = 1.0; //相鄰迭代的結果差n double eps = 1.0e-3; //迭代精度n for(int k=1;(k<Max)&&(dif>eps);k++)n {n dif = 0.0;n printf("n第%d次迭代的結果：n",k);n n for(int i=0;i<size;i++)n {n for(int j=0;j<size;j++)n {n if(i!=j)n {n sum +=A[i][j]*xx[j];n }n }n x[i] = (b[i]-sum)/A[i][i];n sum = 0.0; n } n residual=0.0;n //計算相鄰迭代的結果差n for(int m=0;m<size;m++)n {n dis=fabs(x[m]-xx[m]);n if(dis>residual)n residual=dis;n }n dif=residual;n //列印第k次的結果n int i;n for(i=0;i<size;i++)n {n printf("%12.8f ",x[i]);n xx[i]=x[i];n }n printf("n與上次計算結果的距離(2範數):%12.8f n",dif);n }n printf("n迭代計算的結果為：n");n int k;n for(k=0;k<size;k++)n {n printf("%12.8f ",xx[k]);n }n printf("n");n return 0;n}n

推薦閱讀：