您的当前位置:首页正文

数值分析上机作业(2)

来源:爱站旅游
导读数值分析上机作业(2)
 一、

数值求解如下正方形域上的Poisson方程边值问题

2u2u22f(x,y)1,0x,y1yx  u(0,y)u(1,y)y(1y),0y10x1u(x,0)u(x,1)0,

二、用椭圆型第一边值问题的五点差分格式得到线性方程组为

4ui,jui1,jui1,jui,j1ui,j1h2fij

1i,jNu0,j?,uN1,j?,ui,0?,ui,N1?,0i,jN1

写成矩阵形式Au=f。其中 A11IAIA22I

IANNv1vu2vNb1bf2bN4114Aii114v1(u1,1,u2,1,...,uN,1)T,v2(u1,2,u2,2,...,uN,2)T,......,vN(u1,N,u2,N,...,uN,N)Tb1h2(f1,1,f2,1,...,fN,1)T?,b2h2(f1,2,f2,2,...,fN,2)T?,......,bNh2(f1,N,f2,N,...,fN,N)T?1h,取N9或99,N1fi,j1,i,j1,2,...,N

三 、编写求解线性方程组Au=f的算法程序,用下列方法编程计算,并比较计算速度。 1. 2. 3. 4. 5. 6.

用Jacobi迭代法求解线性方程组Au=f。 用块Jacobi迭代法求解线性方程组Au=f。

用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。 用块SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。 用Gauss-Seidel迭代法求解线性方程组Au=f。 用块Gauss-Seidel迭代法求解线性方程组Au=f。

则h0.1或0.01四、上机报告要求

1.简述方法的基本原理。 2.程序中要加注释。

3.对程序中的主要变量给出说明。 4.附原程序及计算结果。

5.迭代法效率分析:比较迭代法的效率,记录所需的迭代次数,所花费的计算时间。 6.此题可分组做,每组1-3人。需要注明每人完成的任务。

附参考程序:

1、用Gauss-seidel迭代法求解线性方程组Au=f,具体程序如下: function [u,k]=xsgs(n)

f(2:n+1,2:n+1)=(n+1)^(-2)*2; u=zeros(n+2,n+2); e=0.000000001;

for k=1:1000 %迭代求解 er=0;

for j=2:n+1 for i=2:n+1 Ub=u(i,j);

u(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j))/4; er=er+abs(Ub-u(i,j)); %估计误差 end end

if er/n^22、用块Gauss-Seidel迭代法求解线性方程组Au=f,具体程序及结果如下: function[u,k,]=xsbgs(n)

% xsbgs:用块Gauss-seidel迭代法求解线性方程组A*u=f % u: 方程组的解; k: 迭代次数; n: 非边界点数

% a: 方程组系数矩阵Aii的下对角线元素; b: 方程组系数矩阵Aii的主对角线元素; % c: 方程组系数矩阵Aii的上对角线元素;d: 追赶法所求方程的右端向量; % e: 允许误差界; er:迭代误差; f=2*1/(n+1)^2*ones(n+2,n+2);

a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(n+2,n+2);e=0.00001; for k=1:2000 er=0;

for j=2:n+1

Ub=u(:,j);

d(1:n)=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1) ; x=zg(a,b,c,d); % 用追赶法求解

u(2:n+1,j)=x';

er=er+norm(Ub-u(:,j),1); end

if er/n^2附参考程序:解线性方程组Au=f的Gauss-Seidel迭代法。

>> n=9;

b(2:n+1,2:n+1)=0.02; U=zeros(n+2,n+2); e=0.000000001; tic;

for k=1:1000 %迭代求解 er=0;

for j=2:n+1

for i=2:n+1 Ub=U(i,j);

U(i,j)=(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)+b(i,j))/4;

er=er+abs(Ub-U(i,j)); %估计误差 end

end

if er/n^2

因篇幅问题不能全部显示,请点此查看更多更全内容

Top