数值计算方法实验报告(共18页).docx

精选优质文档-倾情为你奉上数值分析实验报告实验一、解线性方程组的直接方法梯形电阻电路问题利用追赶法求解三对角方程组的方法，解决梯形电阻电路问题：电路中的各个电流，须满足下列线性方程组：设，运用追赶法，求各段电路的电流量。问题分析：上述方程组可用矩阵表示为：问题转化为求解，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 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 -2;b=2 5 5 5 5 5 5 5;c=-2 -2 -2 -2 -2 -2 -2 0d=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.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 (error<e)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 = 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=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 (error<10(-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,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 (error<e)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 = gongetidu(A,b,x );运行结果请输入绝对误差限:10(-3)迭代次数:k= 4方程组的解为:x = 1.0000 -2.0000 3.0000 -2.00001.0000结果分析在选择相同的误差容限时，使用以上三种不同的方法来解方程组，对比迭代次数上的不同（1）Jacboi迭代法42次（2）Gauss-Seidel24次（3）共轭梯度法4次，可以看出在收敛速度上共轭梯度法最快，Jacobi迭代法最慢，对比方程组的精确解x=1 -2 3 -2 1，可以看出在收敛效果上共轭梯度法最好！实验三、插值法龙格现象在代数插值中，为了提高插值多项式对函数的逼近层度，常常增加节点个数，级提高多项式的次数，但由于高次多项式插值不收敛，会产生Runge现象，本实验使用函数为y=11+x2，通过对比Lagrange插值法和分段线性插值法可以观察到Runge现场的产生与消失。实验程序函数fxfunction y = fx( x ) y=1./(1+x*x);endlagrange插值方法function = 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(212);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 x<=s(2); l=(x-s(2)./(s(1)-s(2);else l=0;endcase n+1if x>s(n); l=(x-s(n)./(s(n+1)-s(n);else l=0;endotherwiseif x>=s(k-1)&x<=s(k); l=(x-s(k-1)./(s(k)-s(k-1);elseif 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('拟合曲线','实际曲线');gridon;subplot(212);plot(x,e,'b');legend('误差曲线');gridonend区间等分为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函数建立方程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);endfor 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)<ebreak;end k=k+1; h=h/2;endY(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方向网格数目，步进长度为1CX1=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-1for 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 (dif>difmax) difmax=dif;%取误差的最大值endendendendendsubplot(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+X2)的曲线分析Runge现象产生的原因，以及表明了通过增加节点数目来提高函数的逼近程度是不适合的，因此一般情况下不使用多项式插值；（4）使用数值积分中的Romberg型求积公式绘出钟表发条形状的考纽螺线；（5）在专业问题矩形波导中的电位分布问题的分析中，选择了使用差分的方法来代替微分形式的拉普拉斯方程，以及逐次迭代的思想来缩小实际数值计算中的误差，从而绘出了实际的电场分布平面图。这次数值分析实验，我通过编写Matlab程序实现了课堂上所讲过的很多数值计算方法，编写程序的过程就是理解、分析常用数值计算的方法的过程，加深了我对数值计算方法这门课中所讲方法的理解，最重要的是课本中有不少与我们专业相关的问题电阻网络的电流分析、对称点电荷的电场分布、快速傅里叶变换、锁相环路的频率等，通过编程我掌握了使用数值计算的思想来处理与我们专业实际相关问题的方法，体会到了将理论应用于解决实际问题的乐趣。这学期在数值计算方法这门课上收获很多，老师所采用的讲课加实验的教学方式很成功，让我们可以在课堂上学到东西，在课后用Matlab程序实现理论，以及处理和自己专业相关的问题。不足之处，我个人感觉老师可以在我们作报告讲解自己所做的东西时，让其他同学也可以提问充分参与进来。最后，感谢老师这一个学期以来的对我们的认真负责，您上课的认真、讲解的丰富以及对实验的负责深深地感动着每一个学生。