数值计算方法实验报告

数值分析实验报告实验一、解线性方程组的直接方法——梯形电阻电路问题利用追赶法求解三对角方程组的方法，解决梯形电阻电路问题：电路中的各个电流{1i，2i，…，8i}须满足下列线性方程组：RVii22210252321iii0252432iii0252543iii0252654iii0252765iii0252876iii05287ii设V220V，27R，运用追赶法，求各段电路的电流量。问题分析：上述方程组可用矩阵表示为：00000001481.8520000002520000002520000002520000002520000002520000002520000002287654321iiiiiiii问题转化为求解Axb，8阶方阵A满足顺序主子式(1,2...7)0iAi，因此矩阵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: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);fork=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=[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个以回路电流为独立变量的方程，课本上给出的第八个回路电流方程存在问题，正确的应该是78240ii；或者可以根据电路并联分流的知识，同样可以确定78240ii。正确的处理结果应为：x=8.14814.07412.03701.01850.50930.25460.12730.0637实验二、线性方程组的迭代方法——不同迭代法解线性方程组的对比分别用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法；（3）共轭梯度法解线性方程组1234512101234271912314217351732312112435115xxxxx迭代初始向量取(0)00000Tx；(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(errore)break;endendsprintf('迭代次数:\nk=%d',k)sprintf('方程组的解为:\n')xendMatlab命令窗口输入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=42ans=方程组的解为:x=1.0002-2.00022.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');fork=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=[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.00022.9997-2.00000.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:50r=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;endendsprintf('迭代次数:\nk=%d',k)sprintf('方程组的解为:\n')xendMatlab命令窗口输入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.00003.0000-2.00001.0000结果分析在选择相同的误差容限时，使用以上三种不同的方法来解方程组，对比迭代次数上的不同（1）Jacboi迭代法42次（2）Gauss-Seidel24次（3）共轭梯度法4次，可以看出在收敛速度上共轭梯度法最快，Jacobi迭代法最慢，对比方程组的精确解x=[1-23-21]，可以看出在收敛效果上共轭梯度法最好！实验三、插值法——龙格现象在代数插值中，为了提高插值多项式对函数的逼近层度，常常增加节点个数，级提高多项式的次数，但由于高次多项式插值不收敛，会产生Runge现象，本实验使用函数为y=11+x^2，通过对比Lagrange插值法和分段线性插值法可以观察到Runge现场的产生与消失。实验程序函数fxfunctiony=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;forq=1:n+1;l=1;%求插值基函数fork=1:n+1;ifk~=q;l=l.*(x-s(k))./(s(q)-s(k));elsel=l;endendf=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时，如下图所示，黑色曲线为211+x的实际图形，红色曲线即为采用Lagrange插值方法实际描绘出的曲线，从图中可以观察到明显的Runge现象；Runge现象的是因为选择的插值函数并不收敛的缘故。分段线性插值方法function[]=fenduanxianxingchazhi()clearn=input('区间等分份数：\n');s=[-5+10/n*[0:n]];%%%给定的定点，Rf为给定的函数m=0;forx=-5:0.01:5;ff=0;fork=1:n+1;%%%求插值基函数switchkcase1ifx=s(2);l=(x-s(2))./(s(1)-s(2));elsel=0;endcasen+1ifxs(n);l=(x-s(n))./(s(n+1)-s(n));elsel=0;endotherwiseifx=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));elsel=0;endendend%ff=ff+1./(1+25*s(k)*s(k))*l;%%求插值函数值ff=ff+feval(@fx,s(k))*l;%%求插值函数值endm=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份时，如下图所示，黑色曲线为211+x的实际图形，红色曲线即为采用分段线性插值方法实际描绘出的曲线，从图中可以观察到未出现Runge现象，可以观察误差曲线可以看到在将区间等分为20份的时候，和实际曲线的误差限为0.5X10^（-2）；继续提高等份份数，可以继续降低误差！实验四、数值积分——考纽螺线考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为ssdtatsydtatsx020221sin)(21cos)(曲线关于原点对称。取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:1000k=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:kT(k-j+1,j+1)=(4^j*T(k-j+2,j)-T(k-j+1,j))/(4^j-1);endifabs(T(1,k+1)-T(1,k))ebreak;endk=k+1;h=h/2