《计算实习任务1》由会员分享,可在线阅读,更多相关《计算实习任务1(28页珍藏版)》请在金锄头文库上搜索。
1、J d x 2用有限差分法求解正方形域上的 Poisson 方程边值问题+ = f (X, y) = 1,0 x, y 1 Qy2丿u (0, y) = u (1, y) = y (1 - y),0 y 1u (X,0) = u (X,1) = x(1 - X), 0 X 1解:用五点差分格式对椭圆型第一边值问题得到线性方程组为4u u u u u = h2 fi, ji-1, ji+1, ji, j-1i, j+1ij1i, j N在边界上 u= u = j * h *(1 j * h)0, jN+1, ju = u = i * h *(1 i * h),0 i, j maxmax=er;e
2、ndendendU=A;%将矩阵A的值赋予Uif max U,k,er,t=jacobi(9)u =CoIuhltls 1 through 900. 09000. 16000.21000. 24000. 25000. 24000.21000. 16000. 09000. 14160. 18810. 22380. 24600.25350. 24600.22380. 18810. 16000. 18810. 21710. 24110. 25670. 26210.25670. 24110. 21710.21000.22380.24110- 25680.26760.27130. 26760.25680
3、.24110. 24000. 24600. 25670.26760. 27540. 27820.27540.26760. 25670.25000.25350.26210.27130.27820. 28070.27820.27130.26210. 24000. 24600. 25670. 26760. 27540. 27820. 27540. 26760. 25670.21000. 22380.24110. 25680.26760.27130. 26760.25680.24110. 16000. 18810. 21710. 24110. 25670. 26210.25670. 24110. 21
4、710.09000. 14160. 18810.22380. 24600. 25350. 24600.22380. 188100. 09000. 16000.21000. 24000. 25000. 24000.21000. 1600CoIujuiis 10 through 110.090000. 14160. 09000. 18810. 16000.22380. 21000. 24600. 24000.25350. 25000. 24600. 24000.22380. 21000.18810. 16000. 14160. 09000.09000k =336 er =9.5138e-011 t
5、 =0.0113100.20.1结果13360.50.51确定00.40.3Jj&ied on设置Jacobi9.51.3S3e-011迭代诙数计算讨间迭代误差迭代法步长数(2)用块Jacobi迭代法求解线性方程组Au=f。 先把矩阵中的矩阵 Aii 看成 Jacobi 迭代法中的点 对 k=0,1,2 Dx(k+i)二(L + U)x(k)+ b艮卩 x(k+1)二 B x(k) + gJ其中迭代矩阵 B 二 D-1(L + U)二 I - D-1AJ右端向量g二D-ibiiA x(k) ) / Ai ij j ii j=1j知分量形式 X (k+1) = (b -迟 A x(k) 一工 A
6、 x(k) )/A = (b -工i i ij j ij jj =1j =i +1(i=l,2,.n)第二步由于 A 是三对角阵所以对每一个iiX (k+1)i=(b 一乙 A x( k) 一iij j工A x(k)/Aij jiij=1j=i+1=(b一乙A x(k)/Aiij jiij=1j知(i=l,2,.n)用追赶法。追赶法原理:追赶法:AX = b对A进行LU分解如下:b1a2c1b2ur11ur22cn-1bnrn -1unl1nLY 二 dUx 二 Y 计算分解因子阵u二b , r二c (k二1,2,.,n -1)11 k k对 k=2,3, ,nl = a k/, u = b
7、l ck y uk k k kk1 求解 LY=d, y = d11y = d l y (k = 2,3, . .n.),k k k k 1求解 Ux=Y, xnnx = ( y c x ) / u (k = n 1,.,2,1)kk k k+1 k程序:%块Jacobi迭代法function U,k,er,t=bjacobi(n)%变量说明:%u-方程组的解%h-步长%k-迭代次数%n-非边界点数%p-n+2 维向量%q-n+2 维向量%f线性方程组A*u=f的右端矩阵f%a-方程组系数矩阵的下对角线元素%b-方程组系数矩阵的主对角线元素%c-方程组系数矩阵的上对角线元素%d-追赶法所求方程
8、的右端向量%e-允许误差界%er-迭代误差%l -系数矩阵A所分解成的下三角阵L中的下对角线元素l(i)%z -系数矩阵A所分解成的上三角阵U中的主对角线元素z(i)h=1/(n+1);f(2:n+l,2:n+l)= h2;% 初始化 fU=zeros(n+2);%初始化U为n+2阶零矩阵e=0.000000001;%设置误差界a=-1*ones(1,n);b=4*ones(1,n); c=-1*ones(1,n);for j=0:n+1%设置矩阵A, U边界点的值U(1,j+1)=j*h*(1-j*h);U(n+2,j+1)=U(1,j+1);U(j+1,1)=j*h*(1-j*h); U(j+1,n+2)=U(j+1,1);endA=U; tic;% 迭代求解for k=1:1000er=0;for j=2:n+1q=zeros(1,n); q(1)=U(1,j);q(n)=U(n+2,j); d=f(2:n+l,j)+U(2:n+l,j-l)+U(2:n+l,j+l)+q;% 块Jacobi迭代d=d;A(2:n+1,j)=zg(a,b,c