数值计算方法实验报告.docx
《数值计算方法实验报告.docx》由会员分享，可在线阅读，更多相关《数值计算方法实验报告.docx（21页珍藏版）》请在冰豆网上搜索。
数值计算方法实验报告
数值分析实验报告
实验一、解线性方程组的直接方法——梯形电阻电路问题
利用追赶法求解三对角方程组的方法，解决梯形电阻电路问题：
电路中的各个电流{
，
，…，
}须满足下列线性方程组：
设
，
，运用追赶法，求各段电路的电流量。
问题分析：
上述方程组可用矩阵表示为：
问题转化为求解
，8阶方阵A满足顺序主子式
，因此矩阵A存在唯一的Doolittle分解，可以采用解三对角矩阵的追赶法！
追赶法
a=[0-2-2-2-2-2-2-2];
b=[25555555]；
c=[-2-2-2-2-2-2-20];
d=[220/270000000];
Matlab程序
functionx=zhuiganfa（a,b,c,d）
%追赶法实现要求：
|b1|>|C1|>0,|bi|>=|ai|+|ci|
n=length（b）;
u=ones（1,n）;
L=ones（1,n）;
y=ones（1,n）;
u
（1）=b
（1）;
y
（1）=d
（1）;
fori=2:
n
L（i）=a（i）/u（i-1）;
u（i）=b（i）-c（i-1）*L（i）;
y（i）=d（i）-y（i-1）*L（i）;
end
x（n）=y（n）/u（n）;
fork=n-1:
-1:
1
x（k）=（y（k）-c（k）*x（k+1））/u（k）;
end
end
MATLAB命令窗口输入：
a=[0-2-2-2-2-2-2-2];
b=[25555555];
c=[-2-2-2-2-2-2-20]d=[220/270000000];
x=zhuiganfa（a,b,c,d）
运行结果为：
x=
8.14784.07372.03651.01750.50730.25060.11940.0477
存在问题
根据电路分析中的所讲到的回路电流法，可以列出8个以回路电流为独立变量的方程，课本上给出的第八个回路电流方程存在问题，正确的应该是
；或者可以根据电路并联分流的知识，同样可以确定
。
正确的处理结果应为：
x=
8.14814.07412.03701.01850.50930.25460.12730.0637
实验二、线性方程组的迭代方法——不同迭代法解线性方程组的对比
分别用
（1）Jacobi迭代法；
（2）Gauss-Seidel迭代法；（3）共轭梯度法解线性方程组
迭代初始向量取
；
（1）Jacobi迭代法
function[x,k]=jacobi（A,b,x）
%x为迭代初始向量
e=input（'请输入绝对误差限:
\n'）;
D=diag（diag（A））;
n=size（b）;
I=eye（n,n）;
B=I-D\A;
g=D\b;
fork=1:
50;%最大迭代次数为50次
x=B*x+g;
error=norm（b-A*x）;%2-范数
if（errorbreak;
end
end
sprintf（'迭代次数:
\nk=%d',k）
sprintf（'方程组的解为:
\n'）
x
end
Matlab命令窗口输入
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12-2714-1712]';
x=[00000]';
[x,k]=jacobi（A,b,x）;
运行结果
请输入绝对误差限:
10^（-3）
迭代次数:
k=42
ans=
方程组的解为:
x=
1.0002
-2.0002
2.9997
-1.9999
0.9998
（2）Gauss-Seidel迭代法
function[x,k]=gau_sei（A,b,x）
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
e=input（'请输入绝对误差限:
\n'）;
fork=1:
50;
x=（D-L）\U*x+（D-L）\b;
error=norm（b-A*x）;
if（error<10^（-3））
break;
end
end
sprintf（'迭代次数:
\nk=%d',k）
sprintf（'方程组的解为:
\n'）
x
end
Matlab命令窗口输入
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12-2714-1712]';
x=[00000]';
[x,k]=gau_sei（A,b,x）;
运行结果
请输入绝对误差限:
10^（-3）
迭代次数:
k=24
方程组的解为:
x=
1.0002
-2.0002
2.9997
-2.0000
0.9998
（3）共轭梯度法
function[x,k]=gongetidu（A,b,x）
e=input（'请输入绝对误差限:
\n'）;
d=b-A*x;
r=d;
l=（d'*r）./（（A*d）'*d）;
x=x+l*d;
fork=1:
50
r=b-A*x;
t=-（（A*d）'*r）./（（A*d）'*d）;
d=r+t*d;
l=（d'*r）./（（A*d）'*d）;
x=x+l*d;
error=norm（b-A*x）;
if（errorbreak;
end
end
sprintf（'迭代次数:
\nk=%d',k）
sprintf（'方程组的解为:
\n'）
x
end
Matlab命令窗口输入
A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12-2714-1712]';
x=[00000]';
[x,k]=gongetidu（A,b,x）;
运行结果
请输入绝对误差限:
10^（-3）
迭代次数:
k=4
方程组的解为:
x=
1.0000
-2.0000
3.0000
-2.0000
1.0000
结果分析
在选择相同的误差容限时，使用以上三种不同的方法来解方程组，对比迭代次数上的不同
（1）Jacboi迭代法42次
（2）Gauss-Seidel24次（3）共轭梯度法4次，可以看出在收敛速度上共轭梯度法最快，Jacobi迭代法最慢，对比方程组的精确解x=[1-23-21]，可以看出在收敛效果上共轭梯度法最好！
实验三、插值法——龙格现象
在代数插值中，为了提高插值多项式对函数的逼近层度，常常增加节点个数，级提高多项式的次数，但由于高次多项式插值不收敛，会产生Runge现象，本实验使用函数为
，通过对比Lagrange插值法和分段线性插值法可以观察到Runge现场的产生与消失。
实验程序
函数fx
functiony=fx（x）
y=1./（1+x*x）;
end
lagrange插值方法
function[]=lagrangechazhi（）
n=input（'区间等分份数：
\n'）;
s=[-5+10/n*[0:
n]];%%%给定的定点，fx为给定的函数
x=-5:
0.01:
5;
f=0;
forq=1:
n+1;
l=1;%求插值基函数
fork=1:
n+1;
ifk~=q;
l=l.*（x-s（k））./（s（q）-s（k））;
else
l=l;
end
end
f=f+feval（@fx,s（q））*l;%求插值函数
end
y=1./（1+x.*x）;
e=y-f;%误差
subplot（211）;
plot（x,f,'r',x,y,'k'）%作出插值函数曲线
legend（'拟合曲线','实际曲线'）;
gridon;
subplot（212）;
plot（x,e,'b'）;
legend（'误差曲线'）;
gridon
end
区间等分份数为10时，如下图所示，黑色曲线为
的实际图形，红色曲线即为采用Lagrange插值方法实际描绘出的曲线，从图中可以观察到明显的Runge现象；Runge现象的是因为选择的插值函数并不收敛的缘故。
分段线性插值方法
function[]=fenduanxianxingchazhi（）
clear
n=input（'区间等分份数：
\n'）;
s=[-5+10/n*[0:
n]];%%%给定的定点，Rf为给定的函数
m=0;
forx=-5:
0.01:
5;
ff=0;
fork=1:
n+1;%%%求插值基函数
switchk
case1
ifx<=s
（2）;
l=（x-s
（2））./（s
（1）-s
（2））;
else
l=0;
end
casen+1
ifx>s（n）;
l=（x-s（n））./（s（n+1）-s（n））;
else
l=0;
end
otherwise
ifx>=s（k-1）&x<=s（k）;
l=（x-s（k-1））./（s（k）-s（k-1））;
elseifx>=s（k）&x<=s（k+1）;
l=（x-s（k+1））./（s（k）-s（k+1））;
else
l=0;
end
end
end
%ff=ff+1./（1+25*s（k）*s（k））*l;%%求插值函数值
ff=ff+feval（@fx,s（k））*l;%%求插值函数值
end
m=m+1;
f（m）=ff;
end
%%%作出曲线
x=-5:
0.01:
5;
y=1./（1+x.*x）;
e=y-f;%误差
subplot（211）;
plot（x,f,'r',x,y,'k'）%作出插值函数曲线
legend（'拟合曲线','实际曲线'）;
gridon;
subplot（212）;
plot（x,e,'b'）;
legend（'误差曲线'）;
gridon
end
区间等分为20份时，如下图所示，黑色曲线为
的实际图形，红色曲线即为采用分段线性插值方法实际描绘出的曲线，从图中可以观察到未出现Runge现象，可以观察误差曲线可以看到在将区间等分为20份的时候，和实际曲线的误差限为0.5X10^（-2）；继续提高等份份数，可以继续降低误差！
实验四、数值积分——考纽螺线
考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为
曲线关于原点对称。
取a=1，参数s的变化范围[-5，5]，容许误差限分别是10-6和10-10。
选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。
实验程序（程序中编写的积分方法根据数值积分中常见的Romberg求积公式得到）
function[X,Y]=kaoniuluoxian（）
%命令窗口输入[X,Y]=kaoniuluoxian（）
formatlong;
T=zeros（20,20）;
a=input（'积分下限：
\n'）;
b=input（'积分上限：
\n'）;
e=input（'允许误差限：
\n'）;
c=linspace（a,b,1000）;
X=zeros（1,1000）;
Y=zeros（1,1000）;
fx1=inline（'cos（x.^2/2）'）;%考纽螺线方程
fx2=inline（'sin（x.^2/2）'）;%inline函数建立方程
form=1:
1000
k=1;
h=c（m）-a;
T（1,1）=h/2*（fx1（a）+fx1（c（m）））;%计算积分上限为c（m）时fx1的积分值
fori=1:
20%最大迭代次数20次
x=a+h/2:
h:
c（m）-h/2;
T（k+1,1）=T（k,1）/2+h*sum（fx1（x））/2;
forj=1:
k
T（k-j+1,j+1）=（4^j*T（k-j+2,j）-T（k-j+1,j））/（4^j-1）;
end
ifabs（T（1,k+1）-T（1,k））break;
end
k=k+1;
h=h/2;
end
X（m）=T（1,k+1）;
end
form=1:
1000
k=1;
h=c（m）-a;
T（1,1）=h/2*（feval（fx2,a）+feval（fx2,c（m）））;%feval用来求函数fx在给定点处的值
fori=1:
20%计算积分上限为c（m）时fx2的积分值
x=a+h/2:
h:
c（m）-h/2;
T（k+1,1）=T（k,1）/2+h*sum（feval（fx2,x））/2;
forj=1:
k
T（k-j+1,j+1）=（4^j*T（k-j+2,j）-T（k-j+1,j））/（4^j-1）;
end
ifabs（T（1,k+1）-T（1,k））break;
end
k=k+1;
h=h/2;
end
Y（m）=T（1,k+1）;
end
X=-X;
Y=-Y;
plot（X,Y,'k-',-X,-Y,'k-'）;
gridon;
end
命令窗口输入
[X,Y]=kaoniuluoxian（）;
运行结果
积分下限：
0
积分上限：
5
允许误差限：
10^（-6）
绘出图形如下图所示：
容许误差限为10^（-2）时，所绘出图形如下图：
结果分析
对比上述两图，采用Romberg积分方法容许误差限分别为10^（-6）和10^（-2）时，所绘出的图形形状像钟表的发条，即为考纽螺线的形状，但10^（-2）对应的曲线误差偏大！
实验五、同轴矩形金属波导中静电场的分布问题
问题：
同轴矩形金属波导横截面结构如下图所示，其中内金属壁电位为100V，外金属壁为0V，作出其横截面内电位线分布图；
问题背景：
在一个二维矩形区域D内，电位函数
满足二维拉普拉斯方程
在边界L上各点的电位置已给定，即
。
为了采用差分法求解D内的电位函数，作平行于坐标轴的两组平行线
;
将区域D划分为许多正方形网络，网格的边长h称为步长，两组平行线的交点成为网格的节点。
用
表示节点
处的电位值，利用二元函数的Taylor展开式，可将与
相邻的节点上的电位函数值表示为
式中
表示
量级的微量，将以上四式代入拉普拉斯方程可得
这样在
处的电位值
满足的拉普拉斯方程近似的就可以以一个差分方程来代替。
然后以此计算每点的电位，即利用上述近似拉普拉斯方程来计算，用围绕它的四个点的电位平均值作为它的电位值，当所有的点计算完后，用新值代替旧值就完成了一次迭代运算，然后再进行下一次迭代运算，知道计算的新值和旧值之差小于指定的误差容限为止。
就为简单迭代法！
Matlab程序
function[]=dianweifenbu（）
FI=100;
hx=23;hy=23;%外金属槽横截面沿x、y方向网格数目，步进长度为1
CX1=9;
CX2=15;
CY1=9;
CY2=15;%内金属槽占据网格位置:
x输9-15，y轴9-15（刚好在正中间〉
%以上数据可以按要求更改
DX=1+CX2-CX1;
DY=1+CY2-CY1;
v=ones（hy,hx）*50;%每个网格点迭代初值
mesh1=ones（hy,hx）*2;
v（1,:
）=zeros（1,hx）;%上边界初值
mesh1（1,:
）=zeros（1,hx）;
v（hy,:
）=zeros（1,hx）;%下边界初值
mesh1（hy,:
）=zeros（1,hx）;
v（:
1）=zeros（hy,1）;%左边界初值
mesh1（:
1）=zeros（hy,1）;
v（:
hx）=zeros（hy,1）;%右边界初值
mesh1（:
hx）=zeros（hy,1）;
v（CY1:
CY2,CX1:
CX2）=ones（DX,DY）*FI;%内边界初值
mesh1（CY1:
CY2,CX1:
CX2）=ones（DX,DY）;
k=0;%记录迭代次数
difmax=1.0;
while（difmax>1e-6）%误差限
k=k+1;
difmax=0.0;
fori=2:
hy-1
forj=2:
hx-1
m=mesh1（i,j）;
if（m>=2）
oriv=v（i,j）;%原电势值
v（i,j）=1/4*（v（i-1,j）+v（i,j-1）+v（i+1,j）+v（i,j+1））;%拉普拉斯方程差分式
dif=v（i,j）-oriv;%前后两次迭代之差
dif=abs（dif）;
if（dif>difmax）
difmax=dif;%取误差的最大值
end
end
end
end
end
subplot（1,2,1）;
%figure
（1）;
mesh（v）;%画三维曲面图
axis（[-2,hx+3,-2,hy+3,0,100]）;
title（'等势线三维分布图'）;
subplot（1,2,2）;
%figure
（2）;
contour（v,13）;%画等电位线图
holdon;
x=1:
1:
hx;
y=1:
1:
hy;
[xx,yy]=meshgrid（x,y）;%形成栅格
[Gx,Gy]=gradient（v,0.5,0.5）;%计算梯度
quiver（xx,yy,Gx,Gy,'r'）;%根据梯度画箭头
axis（[-2,hx+3,-2,hy+3]）;
plot（[1,1,hx,hx,1],[1,hy,hy,1,1],'k'）;%外框边界线
plot（[CX1,CX1,CX2,CX2,CX1],[CY1,CY2,CY2,CY1,CY1],'k'）;%内框边界线
title（'横截面等势线分布图'）;
text（CX1+0.6,CY1+（CY2-CY1）/2,'U=100V','fontsize',10）;
text（hx/2,hy+1,'U=0V','fontsize',10）;
text（hx/2,0,'U=0V','fontsize',10）;
text（-1.7,hy/2,'U=0V','fontsize',10）;
text（hx+0.4,hy/2,'U=0V','fontsize',10）;
gridon;
holdoff;
end
命令窗口输入
dianweifenbu（）;
运行结果：
从上图可以看出，在靠近内金属波导时，等势线分布更加密集，同时电力线垂直于金属波导边缘，符合实际同轴金属波导中的场强分布规律！
总结
实验报告中涉及到的内容主要包括《数值计算方法》课中所讲解过的
（1）特殊矩阵的三角分解——使用追赶法解决电路中常见的电阻网络问题，通过列出回路电流方程可以发现方程组的系数矩阵恰好为三对角矩阵，既可以使用追赶法解决；
（2）使用迭代法——Jacobi迭代法、Gauss-Seidel迭代法和共轭梯度法来解线性方程组，通过对同一线性方程组设置相同的迭代初始向量和迭代误差，对比计算出在误差允许范围内的近似解所需要的迭代次数，比较了三种方法收敛速度的快慢；（3）使用三角函数插值的方法解决在《数字信号处理》中的快速傅里叶变换问题和使用Lagrange插值、分段线性插值的方法绘出函数1/（1+X^2）的曲线分析Runge现象产生的原因，以及表明了通过增加节点数目来提高函数的逼近程度是不适合的，因此一般情况下不使用多项式插值；（4）使用数值积分中的Romberg型求积公式绘出钟表发条形状的考纽螺线；（5）在专业问题——矩形波导中的电位分布问题的分析中，选择了使用差分的方法来代替微分形式的拉普拉斯方程，以及逐次迭代的思想来缩小实际数值计算中的误差，从而绘出了实际的电场分布平面图。
这次数值分析实验，我通过编写Matlab程序实现了课堂上所讲过的很多数值计算方法，编写程序的过程就是理解、分析常用数值计算的方法的过程，加深了我对《数值计算方法》这门课中所讲方法的理解，最重要的是课本中有不少与我们专业相关的问题——电阻网络的电流分析、对称点电荷的电场分布、快速傅里叶变换、锁相环路的频率等，通过编程我掌握了使用数值计算的思想来处理与我们专业实际相关问题的方法，体会到了将理论应用于解决实际问题的乐趣。
这学期在《数值计算方法》这门课上收获很多，老师所采用的讲课加实验的教学方式很成功，让我们可以在课堂上学到东西，在课后用Matlab程序实现理论，以及处理和自己专业相关的问题。
不足之处，我个人感觉老师可以在我们作报告讲解自己所做的东西时，让其他同学也可以提问充分参与进来。
最后，感谢老师这一个学期以来的对我们的认真负责，您上课的认真、讲解的丰富以及对实验的负责深深地感动着每一个学生。
升级会员 