数值分析上机作业3.doc

1、实验一实验步骤：1、121,1121 21 2222,2 1,1 1,1000nnn n nnnnnnaaaaaA a 得到12naDa21,1,21,100nnnaLa 121,1221,00nnnaUa 则有： ADLUJacobi 迭 代 法 ： ()()()xxbDLUxbDLUxb令 则 称 为 雅 克 比 迭 代 矩 阵)JLfDbJ可 得 雅 克 比 迭 代 的 迭 代 格 式 如 下 ：(0)1(),01,2kkxJf初 始 向 量Gauss-Seidel 迭 代 法 ：()()()()AxbDLUxbDLxUbxDLUxbL令 ， 则 称 为 Gauss-Seidel 迭 代

2、 矩 阵()Gfb可 得 Gauss-Seidel 迭 代 的 迭 代 格 式 如 下 ：(0)1(),01,2kkxf初 始 向 量SOR 迭 代 法 ：0111(1)()()wADLUwLDUwLDwU令 , 则 有 ：()MNAMNxbxbMNxbxbA令 带 入 的 值 可 有WLfb,11(1)()()()wDUwLDUf bL称 为 SOR 迭 代 矩 阵W可 得 SOR 迭 代 的 迭 代 格 式 如 下 ：(0)1(),01,2kkWxLf初 始 向 量2、 编 写 程 序 ：一 、 Jacobi 迭 代 法 M 文 件 ：function y,n=Jacobi(A,b,x0,

3、eps) n=length(A);if nargin=1disp(谱半径大于等于 1，迭代不收敛，无法进行);return;endn=1;for i=1:1:nif sum(A(i,i)=n)=nerror(输入的 A 矩阵的对角线元素不能为 0);return;endendy=B*x0+f;while norm(y-x0)=epsreturn;endn=1;for i=1:1:nif sum(A(i,i)=n)=nerror(输入的 A 矩阵的对角线元素不能为 0);return;endendy=B*x0+f;while norm(y-x0)=eps x0=zeros(20,1); eps=

4、0.0001; y,n=Jacobi(A,b,x0,eps)输出结果：y =0.9637498798116951.1473967056195851.2662264826642621.3048440143034291.3225381929111201.3292674300234981.3320862530763431.3332036239451641.3336501710394001.3338046691840191.3338046691840201.3336501710394001.3332036239451641.3320862530763431.3292674300234981.32253

5、81929111201.3048440143034291.2662264826642621.1473967056195850.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.9637486489082621.1473966333095971.2662281446511131.3048476504060551.3225438351105851.3292750275077231.3320957223

6、809201.3332148363216551.3336629248158521.3338186712192111.3338195394752001.3336654654137911.3332188643199731.3321009563126731.3292811300329621.3225504593923271.3048544570537591.2662347916265421.1474025499265870.963753324289976n =11从结果可以看出其是收敛，在较少的迭代次数下可以满足误差要求的解。（2）第一次给定初始向量 为 20 行一列的 0，右端面项向量 b=20

7、行一列的 1，迭代误差要求)0(x0.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.888477

8、802315454-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(x0.00005，松弛因子为 1.2。命令窗口中输入： b=ones(20,1); x0=zeros(20,1

9、); w=1.2; eps=1e-5; y,n=SOR(A,b,x0,w,eps)输出结果：y =0.2573262365641480.2952446537432390.3215890375263820.3285296272897470.3316773754344230.3327010675976670.3331058406033790.3332480646107180.3333018035055430.3333214781330460.3333272930317400.3333302278521100.3333487204840580.3333089209128300.33317939602

10、79410.3339570317863590.3338959625571420.3243885492191590.3455627826632790.417959490105540n =36从结果，我的得出的结论是松驰系数 在 SOR 迭代中起着相当重要的作用，不同的松驰系数，可能对迭代结果带来很大的影响，不恰当的松驰系数 选取，则可能会得导致无法获得理想的结果，甚至还可能影响到迭代的收敛性。实验二程序设计：编写求解多项式拟合的 Matlab 函数子程序实验要求：用最小二乘法处理下面的实验数据.xi 3 4 5 6 7 8 9fi 2.01 2.98 3.50 5.02 5.47 6.02 7.

11、05并作出 的近似分布图。)(分别采用一次，二次、五次和偶数次多项式来拟合数据得到相应的拟合多项式，并分别作出它们的曲线图。实验步骤：（1）程序：多项式函数 写成 function 的方式，如下jjxfunction y=fai(x,j)y=1;for i=1:jy=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

12、;t=t+fai(s(i),j)*fai(s(i),k);endg(j+1)=t;endA(k+1,:)=g;t=0;for i=1:n;t=t+f(i)*fai(s(i),k);endb(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；五

13、次：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 校正法求解下面非线性方程组1212333cos()0.58sin1.60()xxe（1） 写出 MATLAB 源代码；（2） 给出迭代五次以上的结果；（3） 尝试不同的初值，如可取 ）；(0.1,.)（4） 计算两种方法的用时。实验步

14、骤：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-10x1=x0-A*Fx; k=k+1;if k=5disp(迭代五次的结果为：) ;

15、x1elseif k=100disp(迭代次数过多，不收敛!);return;endFx = subs(F,x,y,z,x1);dFx = subs(dF,x,y,z,x1);A=inv(dFx);x0=x1;enddisp(非线性方程组的解为：)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 = su

16、bs(F,x,y,z,x0);dF = jacobian(F); dFx = subs(dF,x,y,z,x0);A=inv(dFx);k=0;while norm(Fx)1e-10x1=x0-A*Fx;k=k+1;if k=5disp(迭代五次的结果为：) ;x1elseif k=100disp(迭代次数过多，不收敛!);return;endFx1 = 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;enddisp(非线性方程组的解为：)x0 （2）取 ，迭代 20 次的结果：(0.1,.)命令窗口输入：clear,clcx0=0.1,0.1,-0.1;ticmyNewton(x0)t1=tocticmybroyden(x0)t2=toc输出结果：非线性方程组的解为：x0 =0.500000000000000-0.000000000000000-0.523598775598299t1 =0.138421845352998非线性方程组的解为：x0 =0.5000000000003340.000000000000535-0.523598775599102t2 =0.085525802231048