武汉大学数值实习MATLAB作业期末报告.doc

实 验 报 告实验课程：MATLAB数值分析学生姓名： 学 号： 20103010 学 院： 专 业： 指导老师： 2011年12月28日一、问题描述考虑如下Dirichlet问题 其中为正方形区域的边界，类似于模型问题，我们得到差分方程得到系数矩阵A=其中B=-I，I为n-1阶单位矩阵S为对角元为1+次对角元为的n-1阶的对称三对角矩阵。对=sin（xy），=，n=20.用超松弛以及共轭梯度求解差分方程的解。二、Matlab代码及其结果Gauss-Seidel方法：n=20;h=1/n;a=h*h/4;m=n*n;S=zeros(n-1,n-1);for i=1:n-1 S(i,i)=1+h*h/4;endfor i=2:n-1 S(i,i-1)=-1/4;endfor i=1:n-2 S(i,i+1)=-1/4;endB=zeros(n-1,n-1);for i=1:n-1 B(i,i)=-1/4;endA=zeros(n-1)*(n-1),(n-1)*(n-1);for i=1:n-1 A(i-1)*(n-1)+1:i*(n-1),(i-1)*(n-1)+1:i*(n-1)=S;endfor i=1:n-2 A(i*(n-1)+1:(i+1)*(n-1),(i-1)*(n-1)+1:i*(n-1)=B;endfor i=1:n-2 A(i-1)*(n-1)+1:i*(n-1),i*(n-1)+1:(i+1)*(n-1)=B;Endb=1:(n-1)*(n-1)*0;j=1; b(1)=a*sin(1/m)+2*a;b(n-1)=a*sin(n-1)/m)+a*(n-1)*(n-1)+a*(1+n*n);for i=2:n-2 b(i)=a*sin(i*j/m)+a*i*i;endfor j=2:n-2 b(j-1)*(n-1)+1)=a*sin(1*j/m)+a*j*j; b(j*(n-1)=a*sin(n-1)*j/m)+a*(j*j+n*n); for i=(j-1)*(n-1)+2:j*(n-1)-1 b(i)=a*sin(i-(j-1)*(n-1)*j/m); endendj=n-1;b(n-2)*(n-1)+1)=a*sin(1*j/m)+(n-1)*(n-1)*a+(1+n*n)*a;for i=(n-2)*(n-1)+2:(n-1)*(n-1)-1 b(i)=a*sin(i-(n-2)*(n-1)*j/m)+1/4+(i-(n-2)*(n-1)*(n-1)*a;endb(n-1)*(n-1)=a*sin(n-1)*(n-1)/m)+(n*n+(n-1)*(n-1)/(m*2);y=Abn=length (b);Y=zeros (n,1);Y0=Y;Y1=Y;w=1.7 U=zeros (n,n);for i=1:n for j=i+1:n U (i,j)=A (i,j); endendDL=A-U;for II=1:100000 B=-U*Y0+b; Y1 (1)=B (1)/DL(1,1); for i=2:n Y1 (i)=(B (i)-DL (i,1:i-1)*Y1 (1:i-1)/DL(i,i); end Y=w*Y1+(1-w)*Y0; Y0=Y1; r=A*Y-b; err=norm (r); if erreps k=k+1 if k=1 p=r else beta=normr2/normr; p=r+beta*p; normr=normr2; end alpha=(normr)/ (p*A* p); Y=Y+alpha*p; r=r-alpha*A*p; normr2=r*r; k; normr2;endsqrt_normr=sqrt (normr2)krY reshape(Y,19,19) 结果图片： 三、C语言代码及其结果Gauss-Seidel方法：#include#include#includeint main() int i,j,k; int m,n,p; double b400=0,A400400=0,U400=0; double P400=0,r400=0,B400=0; double a; double h; double eps; double normr=0.0,normr2,alpha,beta; double sum,sumr,sumA,sqrt_normr; printf(n= n); scanf(%d,&n); printf(eps= n); scanf(%lf,&eps); printf(%lfn,eps); h=1.0/n; m=n*n; p=n-1; printf(%lfn,h); a=h*h/4; printf(%d,%d,%lf,%lfn,m,p,h,a); j=1; b1=a*sin(1.0/m)+2*a; bn-1=a*sin(1.0*p/m)+a*(p*p)+a*(1+n*n); for(i=2;i=n-2;i+) bi=a*sin(1.0*i*j/m)+a*i*i; for(j=2;j=n-2;j+) b(j-1)*p+1=a*sin(1.0*j/m)+a*j*j; bj*p=a*sin(1.0*p*j/m)+a*(j*j+n*n); for(i=(j-1)*p+2;i=j*p-1;i+) bi=a*sin(i-(j-1)*p)*j*1.0/m); /*for(i=17*p+2;i=18*p;i+) printf(%0.10lfn,bi);*/ j=n-1; b(p-1)*p+1=a*sin(1.0*j/m)+p*p*a+a*(1+n*n); for(i=(p-1)*p+2;i=p*p-1;i+) bi=a*sin(1.0*(i-(p-1)*p)*j/m)+1.0/4+a*(i-(p-1)*p)*(i-(p-1)*p); bp*p=a*sin(1.0*p*p/m)+(n*n+p*p)*2*a; A11=1+a; A12=-0.25; Ap*pp*p-1=-0.25; Ap*pp*p=1+a; for(i=2;i=p*p-1;i+) Aii-1=-0.25; Aii=1+a; Aii+1=-0.25; for(i=n;i=p*p;i+) Ai-n+1i=-0.25; Aii-n+1=-0.25; for(j=2;j=n-1;j+) A(j-1)*p(j-1)*p+1=0.0; A(j-1)*p+1(j-1)*p=0.0; for(i=1;i=p*p;i+) sum=0.0; for(j=1;j=p*p;j+) sum=sum+Aij*Uj; ri=bi-sum; for(i=1;i=p*p;i+) normr=normr+ri*ri; k=0; do k=k+1; if(k=1) for(i=1;i=p*p;i+) Pi=ri; else beta=normr2/normr; for(i=1;i=p*p;i+) Pi=ri+beta*Pi; normr=normr2; sumr=0.0; /r的范数/ for(i=1;i=p*p;i+) sumr=sumr+ri*ri; for(i=1;i=p*p;i+) /求A*P/ Bi=0.0; for(j=1;j=p*p;j+) Bi=Bi+Aij*Pj; sumA=0.0; /A关于P的共轭乘积/ for(i=1;i=p*p;i+) sumA=sumA+Pi*Bi; alpha=sumr/sumA; /求alpha/ for(i=1;i=p*p;i+) Ui=Ui+alpha*Pi; for(i=1;i=p*p;i+) ri=ri-alpha*Bi; normr2=0.0; for(i=1;ieps); printf(U=n); for(i=1;i=p*p;i+) printf(%lfn,Ui); printf(sqrt_normr=%lfn,sqrt_normr); printf(k=%dn,k); system(pause); return 0;SOR方法：#include#include#includeint main() int i,j,k; int n,m,p; double a,h,eps; double w,sum,norm,sumB,sumr; double u400=0,r400=0,B400=0; double u1400=0,u0400=0,b400=0; double A400400=0; printf(n= n); scanf(%d,&n); printf(eps= w= n); scanf(%lf,&eps); scanf(%lf,&w); printf(eps=%lfn,eps); printf(w=%lfn,w); h=1.0/n; m=n*n; p=n-1; printf(%lfn,h); a=h*h/4; printf(%d,%d,%lf,%lfn,m,p,h,a); j=1; b1=a*sin(1.0/m)+2*a; bn-1=a*sin(1.0*p/m)+a*(p*p)+a*(1+n*n); for(i=2;i=n-2;i+) bi=a*sin(1.0*i*j/m)+a*i*i; for(j=2;j=n-2;j+) b(j-1)*p+1=a*sin(1.0*j/m)+a*j*j; bj*p=a*sin(1.0*p*j/m)+a*(j*j+n*n); for(i=(j-1)*p+2;i=j*p-1;i+) bi=a*sin(i-(j-1)*p)*j*1.0/m); j=n-1; b(p-1)*p+1=a*sin(1.0*j/m)+p*p*a+a*(1+n*n); for(i=(p-1)*p+2;i=p*p-1;i+) bi=a*sin(1.0*(i-(p-1)*p)*j/m)+1.0/4+a*(i-(p-1)*p)*(i-(p-1)*p); bp*p=a*sin(1.0*p*p/m)+(n*n+p*p)*2*a; A11=1+a; A12=-0.25; Ap*pp*p-1=-0.25; Ap*pp*p=1+a; for(i=2;i=p*p-1;i+) Aii-1=-0.25; Aii=1+a; Aii+1=-0.25; for(i=n;i=p*p;i+) Ai-n+1i=-0.25; Aii-n+1=-0.25; for(j=2;j=n-1;j+) A(j-1)*p(j-1)*p+1=0.0; A(j-1)*p+1(j-1)*p=0.0; k=0; do k=k+1; for(i=1;i=p*p;i+) /求B=-U*u0+b U用A代替 sumB=0.0; for(j=i+1;j=p*p;j+) sumB=sumB+Aij*u0j; Bi=-sumB+bi; u11=1.0*B1/A11; for(i=2;i=p*p;i+) /求u1 sum=0.0; for(j=1;j=i-1;j+) sum=sum+Aij*u1j; u1i=1.0*(Bi-sum)/Aii; for(i=1;i=p*p;i+) /求解u ui=w*u1i+(1-w)*u0i; for(i=1;i=p*p;i+) /求u0 u0i=u1i; for(i=1;i=p*p;i+) /求残量r sumr=0.0; for(j=