数值分析上机作业第二次.docx

1、数值分析上机作业第二次第二次上机作业数分五班程序分配1-2题：学生姓名：马羽龙 班 级：开发研11-1班学 号：20112120793-4题：学生姓名：马纪翔 班 级：开发研11-1班学 号：20112120965-6题：学生姓名：胡文胜 班 级：化学工艺研11-2班学 号：2011213086一、问题：正方形域上的Poisson方程边值问题，已知条件表述如下：对于上述问题椭圆型第一边值问题的五点差分格式得到线性方程组为：上述线性方程可以写成矩阵形式Au=f。其中对于上述线性方程组可以采用迭代方法求解。二、程序1.用Jacobi迭代法求解线性方程组Au=f。原理：对方程Au=f，Jacobi的

2、迭代式可表达为：程序：function u,k=jacobidd(n)% jacobidd:用jacobi迭代法求解线性方程组A*u=f% u: 方程组的解； k: 迭代次数； n: 非边界点数 % e: 允许误差界； er:迭代误差；%h:迭代步长f(2:n+1,2:n+1)=(n+1)(-2)*2;%根据已知条件，对Au=f中的f矩阵赋值u=zeros(n+2,n+2);%对方程组的解初始化h=1/(n+1); %根据n确定步长hfor j=1:n+2 %利用已知条件计算边界点的值 u(1,j)=(j-1)*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;e

3、nde=0.000000001;%设定允许误差界限值ticfor k=1:1000 %迭代求解 er=0; %初始误差取0 v(1:n+2,1:n+2)= u(1:n+2,1:n+2);%v为中间过渡矩阵 for j=2:n+1 for i=2:n+1 Ub=u(i,j);% 过渡变量，便于计算误差 u(i,j)=(v(i-1,j)+v(i+1,j)+v(i,j-1)+v(i,j+1)+f(i,j)/4;%k次迭代后u(i,j)的值er=er+abs(Ub-u(i,j); %估计误差end end if er/n2 u,k=jacobidd(9)t = 0.0054u = 0 0.0900 0

4、.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.240

5、0 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.090

6、0 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0k = 318 2.用块Jacobi迭代法求解线性方程组Au=f原理：从u的系数矩阵A来看，A为三对角矩阵，可用块迭代法求解Au=f的解，其中，块Jacobi的迭代格式为：对这样格式的方程可采用追赶方法进行求解。程序：function u,k=kjacobidd(n)% kjacobidd:用块jacobi迭代法求解线性方程组A*u=f% u:

7、方程组的解； k: 迭代次数； n: 非边界点数 % a: 方程组系数矩阵 的下对角线元素； b: 方程组系数矩阵 的主对角线元素；% c: 方程组系数矩阵 的上对角线元素；d: 追赶法所求方程的右端向量；% e: 允许误差界； er:迭代误差；%h:迭代步长f=2*1/(n+1)2*ones(n+2,n+2);%根据已知条件，对Au=f中的f矩阵赋值a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(n+2,n+2);%对追赶法中三对角矩阵中分离出来的三个向量赋值e=0.00001;%设定允许误差界限值h=1/(n+1);%根据n确定步

8、长hfor j=1:n+2%利用已知条件计算边界点的值 u(1,j)=(j-1)*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;endfor j=2:n+1 %靠近边界两列单元格的数据 f(2,j)=2*h*h+(j-1)*h*(1-(j-1)*h); f(n+1,j)=2*h*h+(j-1)*h*(1-(j-1)*h);endticfor k=1:2000 er=0;%初始误差取0 v(1:n+2,1:n+2)= u(1:n+2,1:n+2);%v为中间过渡矩阵 for j=2:n+1 Ub=u(:,j); d(1:n)=f(2:n+1,j)+u(2:n+1

9、,j-1)+u(2:n+1,j+1) ; %求取追赶法所需方程的右端向量 x=zg(a,b,c,d); % 用追赶法求解u(2:n+1,j)=x;%u取为x的转置矩阵 er=er+norm(Ub-u(:,j),1); %误差估计 end if er/n2 u,k=kjacobidd(9)t = 0.0181u = 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900

10、 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2099 0.2399 0.2499 0.2399 0.2100 0.1600 0.0900 0 0 0.0900 0.1599 0.2099 0.2399 0.2499 0.2399 0.2100 0.1600 0.0900 0 0 0.0900 0.1599 0.2099 0.2399 0.2499 0.2399 0.2099 0.1600 0.0900 0 0 0.0900 0.1599 0.2099 0.2399 0.2499 0.2

11、399 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2099 0.2399 0.2499 0.2399 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0k = 43

12、3. 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子原理：SOR迭代的分量形式是：而对于Au=f，迭代的分量形式为：其中，w为松弛因子。在确定最佳松弛因子时，以w为变量，取其范围1.11.8,采用0.01的间隔进行编程，通过迭代步数多少，确定最优的w。程序：确定最佳松弛因子function u=wwsordd(n)% wwsordd:最优松弛因子w求解% u: 方程组的解； k: 迭代次数； n: 非边界点数 % e: 允许误差界； er:迭代误差；%w松弛因子for w=1.1:0.01:1.8 %w取值间隔0.01，循环计算最优值f(2:n+1,2:n+1)=(n+1)(-

13、2)*2;u=zeros(n+2,n+2);h=1/(n+1);for j=1:n+2 u(1,j)=(j-1)*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;ende=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)=(w*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j)/4)+(1-w)*u(i,j);er=er+abs(Ub-u(i,j); %估计误差end end if er/n2e,break;end %判

14、断是否达到计算精度，如果达到则退出循环endfprintf(k=%4.0f,w=%3.2fn,k,w)end迭代程序function u,k=sordd(n,w)% sordd:用松弛迭代方法求解线性方程组A*u=f% u: 方程组的解； k: 迭代次数； n: 非边界点数 % e: 允许误差界； er:迭代误差；%w松弛因子；%h:迭代步长f(2:n+1,2:n+1)=(n+1)(-2)*2; %根据已知条件，对Au=f中的f矩阵赋值u=zeros(n+2,n+2); %u初值赋为0h=1/(n+1); %步长值for j=1:n+2 %利用已知条件计算各边界点的值 u(1,j)=(j-1)

15、*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;ende=0.000000001; %定误差界限值ticfor k=1:1000 %迭代求解 er=0; %误差初值定为0 for j=2:n+1 for i=2:n+1 Ub=u(i,j); %计算误差的中间变量 u(i,j)=(w*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j)/4)+(1-w)*u(i,j); %第k步松弛迭代后，u(i,j)的值 er=er+abs(Ub-u(i,j); %估计误差end end if er/n2 u,k=sordd(9,1.55)t

16、 = 2.4944e-004u = 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.160

17、0 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.

18、2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0k = 354. 用块SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子原理：从u的系数矩阵A来看，A为三对角矩阵，可用块迭代法求解Au=f的解，其中，块SOR的迭代格式为：其中，x为采用块迭代得到的解x的转秩矩阵。在确定最佳松弛因子时，采用和第三个编程相同的方法。程序：确定最佳松

19、弛因子functionu=wwksordd(n)% wwksordd:用块sor迭代法求解线性方程组A*u=f最优w求解% u: 方程组的解； k: 迭代次数； n: 非边界点数 % a: 方程组系数矩阵 的下对角线元素； b: 方程组系数矩阵 的主对角线元素；% c: 方程组系数矩阵 的上对角线元素；d: 追赶法所求方程的右端向量；% e: 允许误差界； er:迭代误差；%w松弛因子for w=1.1:0.01:1.8f=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

20、);e=0.00001;h=1/(n+1);for j=1:n+2 u(1,j)=(j-1)*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;endfor j=2:n+1 f(2,j)=2*h*h+(j-1)*h*(1-(j-1)*h); f(n+1,j)=2*h*h+(j-1)*h*(1-(j-1)*h);endfor 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)=w*x

21、+(1-w)*u(2:n+1,j); er=er+norm(Ub-u(:,j),1); end if er/n2e ,break;end % 判断是否达到计算精度，如果达到则退出循环endfprintf(k=%4.0f,w=%3.2fn,k,w)end迭代程序functionu,k=ksordd(n,w)% kgsdd:用块Gauss-seidel迭代法求解线性方程组A*u=f% u: 方程组的解； k: 迭代次数； n: 非边界点数 % a: 方程组系数矩阵 的下对角线元素； b: 方程组系数矩阵 的主对角线元素；% c: 方程组系数矩阵 的上对角线元素；d: 追赶法所求方程的右端向量；%

22、e: 允许误差界； er:迭代误差；%h:迭代步长；%w松弛因子f=2*1/(n+1)2*ones(n+2,n+2); %根据已知条件，对Au=f中的f矩阵赋值a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(n+2,n+2);%对追赶法中三对角矩阵中分离出来的三个向量赋值e=0.00001; %误差界限值h=1/(n+1); %步长值for j=1:n+2 %利用已知条件计算边界点的值 u(1,j)=(j-1)*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;endfor j=2:n+1 %靠近边界两列

23、单元格的数据 f(2,j)=2*h*h+(j-1)*h*(1-(j-1)*h); f(n+1,j)=2*h*h+(j-1)*h*(1-(j-1)*h);endticfor k=1:2000 er=0; %误差初值定为0 for j=2:n+1 Ub=u(:,j); %u为中间过渡矩阵，便于计算误差 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)=w*x+(1-w)*u(2:n+1,j); %SOR方法 er=er+norm(Ub-u(:,j),1);

24、 %误差估计 end if er/n2 u,k=ksordd(9,1.45)t = 0.0126u = 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2

25、100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1

26、600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0k = 155. 用Gauss-Seidel迭代法求解线性方程组Au=f原理：对方程Au=f，Gauss-Seidel的迭代式可表达为：程序：function u,k=gsdd(n)% gsdd:用Gauss-seidel迭

27、代法求解线性方程组A*u=f% u: 方程组的解； k: 迭代次数； n: 非边界点数 % e: 允许误差界； er:迭代误差；h:迭代步长f(2:n+1,2:n+1)=(n+1)(-2)*2;%根据已知条件，对Au=f中的f矩阵赋值u=zeros(n+2,n+2); %对方程组的解初始化h=1/(n+1); %根据n确定步长hfor j=1:n+2 %利用已知条件计算各边界点的值 u(1,j)=(j-1)*h-(j-1)*h)2; u(n+2,j)=(j-1)*h-(j-1)*h)2;ende=0.000000001;%设定允许误差界限值ticfor k=1:1000 %迭代求解 er=0;

28、 %初始误差取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;%k次迭代后u(i,j)的值er=er+abs(Ub-u(i,j); %估计误差end end if er/n2 u,k=gsdd(9)t = 0.0010u = 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.24

29、00 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2500 0.2400 0.2100 0.1600 0.0900 0 0 0.0900 0.1600 0.2100 0.2400 0.2