数值分析上机作业3.doc

实验一实验步骤：1、 得到则有： Jacobi迭代法：令则称为雅克比迭代矩阵可得雅克比迭代的迭代格式如下：Gauss-Seidel迭代法：令，则称为Gauss-Seidel迭代矩阵可得Gauss-Seidel迭代的迭代格式如下：SOR迭代法：令,则有：令带入的值可有称为SOR迭代矩阵可得SOR迭代的迭代格式如下：2、编写程序：一、Jacobi迭代法M文件：function y,n=Jacobi(A,b,x0,eps) n=length(A);if nargin=1 disp(谱半径大于等于1，迭代不收敛，无法进行); return;endn=1;for i=1:1:nif sum(A(i,i)=n)=n error(输入的A矩阵的对角线元素不能为0); return;endendy=B*x0+f;while norm(y-x0)=eps&n100x0=y;y=B*x0+f;n=n+1;end二、Gauss-Seidel迭代法：function y,n=GaussSeidel(A,b,x0,eps)n=length(A);if nargin=1 disp(谱半径大于等于1，迭代不收敛，无法进行); return;endn=1;for i=1:1:nif sum(A(i,i)=n)=n error(输入的A矩阵的对角线元素不能为0); return;endendy=B*x0+f;while norm(y-x0)=eps&n100x0=y;y=B*x0+f;n=n+1;end三、SOR迭代法：function y,n=SOR(A,b,x0,w,eps)n=length(A);if nargin=eps&n b=2.001*ones(20,1); x0=zeros(20,1); eps=0.0001; y,n=Jacobi(A,b,x0,eps)输出结果：y =0.963749879811695 1.147396705619585 1.266226482664262 1.304844014303429 1.322538192911120 1.329267430023498 1.332086253076343 1.333203623945164 1.333650171039400 1.333804669184019 1.333804669184020 1.333650171039400 1.333203623945164 1.332086253076343 1.329267430023498 1.322538192911120 1.304844014303429 1.266226482664262 1.147396705619585 0.963749879811695n =16Gauss-Seidel迭代法求解： b=2.001*ones(20,1); x0=zeros(20,1); eps=0.0001; y,n=GaussSeidel(A,b,x0,eps)输出结果：y =0.963748648908262 1.147396633309597 1.266228144651113 1.304847650406055 1.322543835110585 1.329275027507723 1.332095722380920 1.333214836321655 1.333662924815852 1.333818671219211 1.333819539475200 1.333665465413791 1.333218864319973 1.332100956312673 1.329281130032962 1.322550459392327 1.304854457053759 1.266234791626542 1.147402549926587 0.963753324289976n =11 从结果可以看出其是收敛，在较少的迭代次数下可以满足误差要求的解。（2） 第一次给定初始向量为20行一列的0，右端面项向量b=20行一列的1，迭代误差要求0.00005，松弛因子为 1.5。 命令窗口输入： b=ones(20,1); x0=zeros(20,1); w=1.5; eps=1e-5; y,n=SOR(A,b,x0,w,eps)输出结果：y =1.0e+012 * -0.508184555095282 -0.968955023439180 -1.540013778883930 -2.173751238422641 -2.876743633728070 -3.635630337032175 -4.437482590297024 -5.263455839480649 -6.090077734793300 -6.888477802315454 -7.624337014885191 -8.257825755757494 -8.743723731466995 -9.031937365363348 -9.067525660623353 -8.793955338248141 -8.145209778158643 -7.083095899003873 -5.459752204553985 -3.565052640053931n =100第二次给定初始向量为20行一列的0，右端面项向量b=20行一列的1，迭代误差要求0.00005，松弛因子为 1.2。命令窗口中输入： b=ones(20,1); x0=zeros(20,1); w=1.2; eps=1e-5; y,n=SOR(A,b,x0,w,eps)输出结果：y =0.257326236564148 0.295244653743239 0.321589037526382 0.328529627289747 0.331677375434423 0.332701067597667 0.333105840603379 0.333248064610718 0.333301803505543 0.333321478133046 0.333327293031740 0.333330227852110 0.333348720484058 0.333308920912830 0.333179396027941 0.333957031786359 0.333895962557142 0.324388549219159 0.345562782663279 0.417959490105540n =36从结果，我的得出的结论是松驰系数在SOR迭代中起着相当重要的作用，不同的松驰系数，可能对迭代结果带来很大的影响，不恰当的松驰系数选取，则可能会得导致无法获得理想的结果，甚至还可能影响到迭代的收敛性。实验二程序设计：编写求解多项式拟合的Matlab函数子程序实验要求：用最小二乘法处理下面的实验数据.xi3456789fi2.012.983.505.025.476.027.05并作出的近似分布图。分别采用一次，二次、五次和偶数次多项式来拟合数据得到相应的拟合多项式，并分别作出它们的曲线图。实验步骤：（1） 程序：多项式函数写成function的方式，如下function y=fai(x,j)y=1;for i=1:j y=x.*y;end主程序：clears=3 4 5 6 7 8 9;f=2.01 2.98 3.50 5.02 5.47 6.02 7.05;%计算给定的数据点的数目n=length(f);%给定需要拟合的数据的最高次多项式的次数m=1;for k=0:m; g=zeros(1,m+1); for j=0:m; t=0; for i=1:n; t=t+fai(s(i),j)*fai(s(i),k); end g(j+1)=t; end A(k+1,:)=g; t=0; for i=1:n; t=t+f(i)*fai(s(i),k); end b(k+1,1)=t;enda=Ab%求出多项式系数x=s(1):0.01:s(n);y=0;for i=0:m; y=y+a(i+1)*fai(x,i);endplot(x,y)grid onhold on plot(s,f,rx) 一次，二次、五次和偶数次多项式来拟合数据得到相应的拟合多项式：一次：y1=-0.38643+0.82750x；二次：y2=-1.03024+1.06893x-0.02012x2；五次：y5=-50.75309+51.53527x-19.65947x2+3.66585x3-0.32886x4+0.01137x5；六次：y6= -514.9848+ 587.1498x -268.9138x2+63.6654x3 -8.2239x4+ 0.5509x5 -0.0150x6；实验三实验内容与要求：分别用Newton法用Broyden秩1校正法求解下面非线性方程组 （1） 写出MATLAB源代码；（2） 给出迭代五次以上的结果；（3） 尝试不同的初值，如可取）；（4） 计算两种方法的用时。实验步骤：Newton法源代码：functiont=myNewton(x0) syms x y z;x0 = transpose(x0);f1=3*x-cos(y*z)-0.5;f2=x2-81*(y+0.1)2+sin(z)+1.06; f3=exp(-x*y)+20*z+1/3*(10*pi-3);F=f1;f2;f3;Fx = subs(F,x,y,z,x0);dF = jacobian(F); dFx = subs(dF,x,y,z,x0);A=inv(dFx);k=0;while norm(Fx)1e-10 x1=x0-A*Fx; k=k+1; if k=5 disp(迭代五次的结果为：) ; x1 elseif k=100 disp(迭代次数过多，不收敛!); return; end Fx = subs(F,x,y,z,x1); dFx = subs(dF,x,y,z,x1); A=inv(dFx); x0=x1;end disp(非线性方程组的解为：) x0 Broyden秩1校正法源代码：functiont=mybroyden(x0)syms x y z;x0 = transpose(x0);f1=3*x-cos(y*z)-0.5;f2=x2-81*(y+0.1)2+sin(z)+1.06; f3=exp(-x*y)+20*z+1/3*(10*pi-3);F=f1;f2;f3;Fx = subs(F,x,y,z,x0);dF = jacobian(F); dFx = subs(dF,x,y,z,x0);A=inv(dFx);k=0;while norm(Fx)1e-10 x1=x0-A*Fx; k=k+1; if k=5 disp(迭代五次的结果为：) ; x1 elseif k=100 disp(迭代次数过多，不收敛!); return; end Fx1 = subs(F,x,y,z,x1); s=x1-x0; r=Fx1-Fx; A=A-(A*r-s)*s*A)/(s*A*r); x0=x1; Fx=Fx1;end disp(非线性方程组的解为：) x0 （2） 取，迭代20次的结果：命令窗口输入：clear,clcx