数值计算方法实验报告.docx

1、数值计算方法实验报告数值分析实验报告实验一、解线性方程组的直接方法梯形电阻电路问题 利用追赶法求解三对角方程组的方法，解决梯形电阻电路问题：电路中的各个电流，须满足下列线性方程组：设，运用追赶法，求各段电路的电流量。问题分析：上述方程组可用矩阵表示为： 问题转化为求解，8阶方阵A满足顺序主子式，因此矩阵A存在唯一的Doolittle分解，可以采用解三对角矩阵的追赶法！追赶法 a=0 -2 -2 -2 -2 -2 -2 -2; b=2 5 5 5 5 5 5 5； c=-2 -2 -2 -2 -2 -2 -2 0; d=220/27 0 0 0 0 0 0 0;Matlab程序 function

2、 x= 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:nL(i)=a(i)/u(i-1);u(i)=b(i)-c(i-1)*L(i);y(i)=d(i)-y(i-1)*L(i);endx(n)=y(n)/u(n);for k=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1)/u(k);endendMATLAB命令窗口输入： a=0 -2 -2 -2 -2 -2 -2

3、-2; b=2 5 5 5 5 5 5 5; c=-2 -2 -2 -2 -2 -2 -2 0 d=220/27 0 0 0 0 0 0 0; x= zhuiganfa(a,b,c,d )运行结果为：x =8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477存在问题 根据电路分析中的所讲到的回路电流法，可以列出8个以回路电流为独立变量的方程，课本上给出的第八个回路电流方程存在问题，正确的应该是；或者可以根据电路并联分流的知识，同样可以确定。正确的处理结果应为：x =8.1481 4.0741 2.0370 1.0185 0.5093 0

4、.2546 0.1273 0.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-DA;g=Db;for k=1:50;%最大迭代次数为50次 x=B*x+g; error=norm(b-A*x);%2-范数if (errore)break;enden

5、dsprintf(迭代次数:nk= %d,k)sprintf(方程组的解为:n )xendMatlab命令窗口输入A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15;b=12 -27 14 -17 12;x=0 0 0 0 0;x,k = jacobi(A,b,x);运行结果请输入绝对误差限:10(-3)迭代次数:k= 42ans =方程组的解为:x = 1.0002 -2.0002 2.9997 -1.99990.9998(2)Gauss-Seidel迭代法function x,k = gau_sei(A,b,x)D=

6、diag(diag(A);L=D-tril(A);U=D-triu(A);e=input(请输入绝对误差限:n);for k=1:50; x=(D-L)U*x+(D-L)b;error=norm(b-A*x);if (error10(-3)break;endendsprintf(迭代次数:nk= %d,k)sprintf(方程组的解为:n )xendMatlab命令窗口输入A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15;b=12 -27 14 -17 12;x=0 0 0 0 0;x,k = gau_sei(A,b,

7、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;for k=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 (errore)break;end

8、endsprintf(迭代次数:nk= %d,k)sprintf(方程组的解为:n )xendMatlab命令窗口输入A=10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15;b=12 -27 14 -17 12;x=0 0 0 0 0;x,k = gongetidu(A,b,x );运行结果请输入绝对误差限:10(-3)迭代次数:k= 4方程组的解为:x = 1.0000 -2.0000 3.0000 -2.00001.0000结果分析 在选择相同的误差容限时，使用以上三种不同的方法来解方程组，对比迭代次数上的不同（1）Ja

9、cboi迭代法42次（2）Gauss-Seidel24次（3）共轭梯度法4次，可以看出在收敛速度上共轭梯度法最快，Jacobi迭代法最慢，对比方程组的精确解x=1 -2 3 -2 1，可以看出在收敛效果上共轭梯度法最好！实验三、插值法龙格现象 在代数插值中，为了提高插值多项式对函数的逼近层度，常常增加节点个数，级提高多项式的次数，但由于高次多项式插值不收敛，会产生Runge现象，本实验使用函数为，通过对比Lagrange插值法和分段线性插值法可以观察到Runge现场的产生与消失。实验程序函数fxfunction y = fx( x ) y=1./(1+x*x);endlagrange插值方法f

10、unction = lagrangechazhi( )n=input(区间等分份数：n);s=-5+10/n*0:n;%给定的定点，fx为给定的函数x=-5:0.01:5;f=0;for q=1:n+1; l=1;%求插值基函数for k=1:n+1;if k=q; l=l.*(x-s(k)./(s(q)-s(k);else l=l;endend f=f+feval(fx,s(q)*l;%求插值函数endy=1./(1+x.*x);e=y-f;%误差subplot(211);plot(x,f,r,x,y,k)%作出插值函数曲线legend(拟合曲线,实际曲线);gridon;subplot(2

11、12);plot(x,e,b);legend(误差曲线);gridonend区间等分份数为10时，如下图所示，黑色曲线为的实际图形，红色曲线即为采用Lagrange插值方法实际描绘出的曲线，从图中可以观察到明显的Runge现象；Runge现象的是因为选择的插值函数并不收敛的缘故。分段线性插值方法function = fenduanxianxingchazhi( )clearn=input(区间等分份数：n);s=-5+10/n*0:n;%给定的定点，Rf为给定的函数m=0;for x=-5:0.01:5;ff=0;for k=1:n+1;%求插值基函数switch kcase 1if xs(n

12、); l=(x-s(n)./(s(n+1)-s(n);else l=0;endotherwiseif x=s(k-1)&x=s(k)&x=s(k+1); l=(x-s(k+1)./(s(k)-s(k+1);else l=0;endendend%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(拟合曲线,实际曲

13、线);gridon;subplot(212);plot(x,e,b);legend(误差曲线);gridonend 区间等分为20份时，如下图所示，黑色曲线为的实际图形，红色曲线即为采用分段线性插值方法实际描绘出的曲线，从图中可以观察到未出现Runge现象，可以观察误差曲线可以看到在将区间等分为20份的时候，和实际曲线的误差限为0.5X10（-2）；继续提高等份份数，可以继续降低误差！实验四、数值积分考纽螺线 考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为曲线关于原点对称。取a=1，参数s的变化范围-5，5，容许误差限分别是10-6和10-10。选取适当的节点个数，利用

14、数值积分方法计算曲线上点的坐标，并画出曲线的图形。 实验程序（程序中编写的积分方法根据数值积分中常见的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=inlin

15、e(sin(x.2/2);%inline函数建立方程for m=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; for j=1:k T(k-j+1,j+1)=(4j*T(k-j+2,j)-T(k-j+1,j)/(4j-1);endif abs(T(1,k+1)-T(1,k)ebreak;end k=k+1; h=h/2;endX(m)=T(1,k+1);en

16、dfor m=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; for j=1:k T(k-j+1,j+1)=(4j*T(k-j+2,j)-T(k-j+1,j)/(4j-1);endif abs(T(1,k+1)-T(1,k)1e-6)%误差限 k=k+1;difmax=0.0;fori=2:hy-1fo

17、r j=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 (difdifmax) difmax=dif;%取误差的最大值endendendendendsubplot(1,2,1);%figure(1);mesh(v);%画三维曲面图axis(-2,hx+3,-2,hy+3,0,100);title(等势线三维分布图);subplot(1,2,2);%figu

18、re(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=

19、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）特殊矩阵的三角分解

20、使用追赶法解决电路中常见的电阻网络问题，通过列出回路电流方程可以发现方程组的系数矩阵恰好为三对角矩阵，既可以使用追赶法解决；（2）使用迭代法Jacobi迭代法、Gauss-Seidel迭代法和共轭梯度法来解线性方程组，通过对同一线性方程组设置相同的迭代初始向量和迭代误差，对比计算出在误差允许范围内的近似解所需要的迭代次数，比较了三种方法收敛速度的快慢；（3）使用三角函数插值的方法解决在数字信号处理中的快速傅里叶变换问题和使用Lagrange插值、分段线性插值的方法绘出函数1/(1+X2)的曲线分析Runge现象产生的原因，以及表明了通过增加节点数目来提高函数的逼近程度是不适合的，因此一般情况下

21、不使用多项式插值；（4）使用数值积分中的Romberg型求积公式绘出钟表发条形状的考纽螺线；（5）在专业问题矩形波导中的电位分布问题的分析中，选择了使用差分的方法来代替微分形式的拉普拉斯方程，以及逐次迭代的思想来缩小实际数值计算中的误差，从而绘出了实际的电场分布平面图。 这次数值分析实验，我通过编写Matlab程序实现了课堂上所讲过的很多数值计算方法，编写程序的过程就是理解、分析常用数值计算的方法的过程，加深了我对数值计算方法这门课中所讲方法的理解，最重要的是课本中有不少与我们专业相关的问题电阻网络的电流分析、对称点电荷的电场分布、快速傅里叶变换、锁相环路的频率等，通过编程我掌握了使用数值计算的思想来处理与我们专业实际相关问题的方法，体会到了将理论应用于解决实际问题的乐趣。 这学期在数值计算方法这门课上收获很多，老师所采用的讲课加实验的教学方式很成功，让我们可以在课堂上学到东西，在课后用Matlab程序实现理论，以及处理和自己专业相关的问题。不足之处，我个人感觉老师可以在我们作报告讲解自己所做的东西时，让其他同学也可以提问充分参与进来。 最后，感谢老师这一个学期以来的对我们的认真负责，您上课的认真、讲解的丰富以及对实验的负责深深地感动着每一个学生。