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

资源描述

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

《数值分析上机作业第二次.docx》由会员分享，可在线阅读，更多相关《数值分析上机作业第二次.docx（17页珍藏版）》请在冰豆网上搜索。

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

数值分析上机作业第二次

第二次上机作业

数分五班

程序分配

1-2题：

学生姓名：

马羽龙

班级：

开发研11-1班

学号：

2011212079

3-4题：

学生姓名：

马纪翔

班级：

开发研11-1班

学号：

2011212096

5-6题：

学生姓名：

胡文胜

班级：

化学工艺研11-2班

学号：

2011213086

一、问题：

正方形域上的Poisson方程边值问题，已知条件表述如下：

对于上述问题椭圆型第一边值问题的五点差分格式得到线性方程组为：

上述线性方程可以写成矩阵形式Au=f。

其中

对于上述线性方程组可以采用迭代方法求解。

二、程序

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

原理：

对方程Au=f，Jacobi的迭代式可表达为：

程序：

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确定步长h

forj=1:

n+2%利用已知条件计算边界点的值

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

e=0.000000001;%设定允许误差界限值

tic

fork=1:

1000%迭代求解

er=0;%初始误差取0

v（1:

n+2,1:

n+2）=u（1:

n+2,1:

n+2）;%v为中间过渡矩阵

forj=2:

n+1

fori=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

ifer/n^2

end

toc

结果：

>>[u,k]=jacobidd（9）

0.0054

00.09000.16000.21000.24000.25000.24000.21000.16000.09000

318

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

原理：

从u的系数矩阵A来看，A为三对角矩阵，可用块迭代法求解Au=f的解，其中，块Jacobi的迭代格式为：

对这样格式的方程可采用追赶方法进行求解。

程序：

function[u,k]=kjacobidd（n）

%kjacobidd:

用块jacobi迭代法求解线性方程组A*u=f

%u:

方程组的解；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确定步长h

forj=1:

n+2%利用已知条件计算边界点的值

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

forj=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）;

end

tic

fork=1:

2000

er=0;%初始误差取0

v（1:

n+2,1:

n+2）=u（1:

n+2,1:

n+2）;%v为中间过渡矩阵

forj=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';%u取为x的转置矩阵

er=er+norm（Ub-u（:

j）,1）;%误差估计

end

ifer/n^2

end

t=toc

结果：

>>[u,k]=kjacobidd（9）

0.0181

00.09000.16000.21000.24000.25000.24000.21000.16000.09000

00.09000.16000.20990.23990.24990.23990.21000.16000.09000

00.09000.15990.20990.23990.24990.23990.21000.16000.09000

00.09000.15990.20990.23990.24990.23990.20990.16000.09000

00.09000.15990.20990.23990.24990.23990.21000.16000.09000

00.09000.16000.20990.23990.24990.23990.21000.16000.09000

00.09000.16000.21000.24000.25000.24000.21000.16000.09000

3.用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子

原理：

SOR迭代的分量形式是：

而对于Au=f，迭代的分量形式为：

其中，w为松弛因子。

在确定最佳松弛因子时，以w为变量，取其范围1.1~1.8,采用0.01的间隔进行编程，通过迭代步数多少，确定最优的w。

程序：

确定最佳松弛因子

function[u]=wwsordd（n）

%wwsordd:

最优松弛因子w求解

%u:

方程组的解；k:

迭代次数；n:

非边界点数

%e:

允许误差界；er:

迭代误差；%w松弛因子

forw=1.1:

0.01:

1.8%w取值间隔0.01，循环计算最优值

f（2:

n+1,2:

n+1）=（n+1）^（-2）*2;

u=zeros（n+2,n+2）;

h=1/（n+1）;

forj=1:

n+2

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

e=0.000000001;

fork=1:

1000%迭代求解

er=0;

forj=2:

n+1

fori=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

ifer/n^2

end

fprintf（'k=%4.0f,w=%3.2f\n',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初值赋为0

h=1/（n+1）;%步长值

forj=1:

n+2%利用已知条件计算各边界点的值

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

e=0.000000001;%定误差界限值

tic

fork=1:

1000%迭代求解

er=0;%误差初值定为0

forj=2:

n+1

fori=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

ifer/n^2

end

t=toc

结果：

由程序得最佳松弛因子为1.55，代入迭代程序有：

>>[u,k]=sordd（9,1.55）

2.4944e-004

00.09000.16000.21000.24000.25000.24000.21000.16000.09000

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

原理：

从u的系数矩阵A来看，A为三对角矩阵，可用块迭代法求解Au=f的解，其中，块SOR的迭代格式为：

其中，x’为采用块迭代得到的解x的转秩矩阵。

在确定最佳松弛因子时，采用和第三个编程相同的方法。

程序：

确定最佳松弛因子

function[u]=wwksordd（n）

%wwksordd:

用块sor迭代法求解线性方程组A*u=f最优w求解

%u:

方程组的解；k:

迭代次数；n:

非边界点数

%a:

方程组系数矩阵的下对角线元素；b:

方程组系数矩阵的主对角线元素；

%c:

方程组系数矩阵的上对角线元素；d:

追赶法所求方程的右端向量；

%e:

允许误差界；er:

迭代误差；%w松弛因子

forw=1.1:

0.01:

1.8

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;

h=1/（n+1）;

forj=1:

n+2

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

forj=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）;

end

fork=1:

2000

er=0;

forj=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'+（1-w）*u（2:

n+1,j）;

er=er+norm（Ub-u（:

j）,1）;

end

ifer/n^2

end

fprintf（'k=%4.0f,w=%3.2f\n',k,w）

end

迭代程序

function[u,k]=ksordd（n,w）

%kgsdd:

用块Gauss-seidel迭代法求解线性方程组A*u=f

%u:

方程组的解；k:

迭代次数；n:

非边界点数

%a:

方程组系数矩阵的下对角线元素；b:

方程组系数矩阵的主对角线元素；

%c:

方程组系数矩阵的上对角线元素；d:

追赶法所求方程的右端向量；

%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）;%步长值

forj=1:

n+2%利用已知条件计算边界点的值

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

forj=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）;

end

tic

fork=1:

2000

er=0;%误差初值定为0

forj=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）;%误差估计

end

ifer/n^2

end

t=toc

结果：

用程序得最优松弛因子w为1.45，带入迭代程序有：

>>[u,k]=ksordd（9,1.45）

0.0126

00.09000.16000.21000.24000.25000.24000.21000.16000.09000

5.用Gauss-Seidel迭代法求解线性方程组Au=f

原理：

对方程Au=f，Gauss-Seidel的迭代式可表达为：

程序：

function[u,k]=gsdd（n）

%gsdd:

用Gauss-seidel迭代法求解线性方程组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确定步长h

forj=1:

n+2%利用已知条件计算各边界点的值

u（1,j）=（j-1）*h-（（j-1）*h）^2;

u（n+2,j）=（j-1）*h-（（j-1）*h）^2;

end

e=0.000000001;%设定允许误差界限值

tic

fork=1:

1000%迭代求解

er=0;%初始误差取0

forj=2:

n+1

fori=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

ifer/n^2

end

t=toc

结果：

>>[u,k]=gsdd（9）

0.0010

00.09000.16000.21000.24000.25000.24000.21000.16000.09000

00.09000.16000.21000.24000.2

展开阅读全文