下载文档原格式
一．实验任务用MA TLAB 语言编写连续函数最佳平方逼近的算法程序（函数式M 文件）。
并用此程序进行数值试验，写出实验报告。
二．实验方法最佳平方逼近方法采用基于正交多项式的最佳平方逼近，选择Lengendre 多项式做基。
计算组合系数时，函数的积分采用变步长复化梯形求积法。
三．程序功能和使用说明1.采用基于正交多项式的最佳平方逼近，选择Lengendre 多项式做基利用递推关系0112()1,()()(21)()(1)()/2,3,.....n n n P x P x xP x n xP x n P x n n --===---⎡⎤⎣⎦=可构造出用户需要的任意次数的最佳平方逼近多项式。
2. 用M 文件建立数学函数，实现程序通过修改建立数学函数的M 文件以适用不同的被逼近函数。
3.已经考虑一般的情况]1,1[],[)(+-≠∈b a x f ，程序有变量代换的功能。
4.计算组合系数时，函数的积分采用变步长复化梯形求积法5.可根据需要，求出二次、三次、。
最佳平方逼近函数）x s (。
6.最后作出逼近函数）x s (和被逼近函数)(x f 的曲线图可进行比较，分别用绘图函数plot 和fplot 绘图。
7.在matlab 的命令窗口，输入[c,sx]=leastp(@func1,a,b,n)，func1是被逼近函数，b 和a 分别是逼近函数的上、下区间，n 为最佳平方逼近的次数，可为任意次数。
四．程序代码（含注释）1. 最佳平方逼近主函数function [c,sx]=leastp(func,a,b,n)%LEASTP.m:least-square fitting with legendre polynomials%func 指被逼近函数，调用需要用句柄%a,b 分别指被逼近函数的区间上下限%n 指最佳平方逼近的次数syms t;syms x;%以Lengendre 多项式为基，构造任意次数的最佳平方逼近多项式p(2)=t;p(1)=1;if n>1for j=3:1:(n+1)p(j)=((2*j-3)*t*p(j-1)-(j-2)*p(j-2))/(j-1);endend%变量代换，区间调整为[-1,1]f=feval(func,(b-a)/2*t+(b+a)/2);%计算组合系数，其中调用变步长复化梯形求积函数trapzfor j=1:1:(n+1)c(j)=(2*j-1)/2*trapz(f*p(j),-1,1);end%将组合系数与对应的最佳平方多项式相乘然后求和，得到最佳逼近函数sx=0;for j=1:1:(n+1)sx=sx+c(j)*p(j);end%将变量替换还原sx=subs(sx,(2*x-a-b)/(b-a));%使用fplot绘制原函数图像f1=feval(func,x);f1=inline(f1);[x,y]=fplot(f1,[a,b]);plot(x,y,'r-','linewidth',1.5);hold on;%使用plot绘制最佳平方逼近函数图像g=linspace(a,b,(b-a)*300);fsx=subs(sx,g);plot(g,fsx,'b-','linewidth',1.5);str=strcat(num2str(n),'次最佳平方逼近');legend('原函数',str);end2. 计算组合系数，变步长复化梯形求积法function To1=trapz(func,a,b)%半分区间复化梯形公式计算定积分%func指需要求积分的原函数%a,b分别指积分上下区间%初值h=b-a;To=(subs(func,a)+subs(func,b))*(b-a)/2;e=1;while e>10^-6%迭代终止条件，前后两次积分值差小于10^-6 H=0;x=a+h/2;while x<bH=H+subs(func,x);%计算出所有二分新出现的值的和x=x+h;endTo1=0.5*(To+h*H);%计算出新的积分值e=abs(To1-To);h=h/2;%继续半分区间，进行迭代计算To=To1;endend3. 以.m文件定义被逼近函数function y=func1(x)y=x*cos(x);end五．实验结果1. 一次最佳平方逼近c =-1.1702 -2.4235sx=1.253290 - 1.211752*x2. 二次最佳平方逼近c =-1.1702 -2.4235 -0.4265sx=-0.159939*x^2 - 0.571997*x + 0.8267873. 三次最佳平方逼近c =-1.1702 -2.4235 -0.4265 1.2216sx=0.381759*x^3 - 2.450495*x^2 + 3.092892*x - 0.3948434. 四次最佳平方逼近c =-1.1702 -2.4235 -0.4265 1.2216 0.3123sx =0.085392*x^4 - 0.301375*x^3 - 0.693864*x^2 + 1.531443*x - 0.082553六．分析与讨论从次数从1到4的最佳平方逼近图像对比可以发现，次数越高，图像拟合效果越好。
数值分析实验报告一、 实验3.1 题目：考虑线性方程组b Ax =，n n R A ⨯∈，n R b ∈，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss 消去过程。
（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=6816816816 A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157 b ，则方程有解()T x 1,,1,1*⋯=。
取10=n 计算矩阵的条件数。
分别用顺序Gauss 消元、列主元Gauss 消元和完全选主元Gauss 消元方法求解，结果如何？（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。
1. 算法介绍首先，分析各种算法消去过程的计算公式， 顺序高斯消去法：第k 步消去中，设增广矩阵B 中的元素()0k kk a ≠（若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k 行以下各行计算()(),1,2,,k ikik k kka l i k k n a ==++，分别用ik l -乘以增广矩阵B 的第k 行并加到第1,2,,k k n ++行，则可将增广矩阵B 中第k 列中()k kka 以下的元素消为零；重复此方法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即()()(),n n n B Ab ⎡⎤=⎣⎦； 列主元高斯消去法：第k 步消去中，在增广矩阵B 中的子方阵()()()()k kkkknk k nknn a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦中，选取()k k i k a 使得()(k)max k ki k ik k i na a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换，然后按照和顺序消去法相同的步骤进行。
《数值分析》计算实习报告第二题院系：机械工程及自动化学院学号：姓名：2017年11月一、题目要求试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知a ij ={sin (0.5i +0.2j ) i ≠j1.52cos (i +1.2j ) i =j(i,j =1,2, (10)说明：1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10−12。
2.打印以下内容： （1）全部源程序；（2）矩阵A 经过拟上三角化后所得的矩阵A (n−1)； （3）对矩阵A (n−1)实行QR 方法迭代结束后所得的矩阵；（4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,⋯,10)，其中R i =Re(λi )，I i = Im(λi ) 。
若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。
3.采用e 型数输出实型数，并且至少显示12位有效数字。
二、算法设计思路和方案1. 将矩阵A 拟上三角化得到矩阵A (n−1)为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n−1)，然后用QR 方法计算A (n−1)的全部特征值，而A (n−1)的特征值就是A 的特征值。
具体算法如下：记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2,,;,1,,)ij a i n j r r n ==+。
对于1,2,,2r n =-执行（1）若()(2,3,,)r ira i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。
（2）计算r d =()()1,sgn r r r r r c a d +=- (若()1,0r r r a +=，则取r r c d =)2()1,r r r r r r h c c a +=-（3）令()()()()1,2,0,,0,,,,T r r r n r r rr r rnru ac aaR ++=-∈。
第二次计算实验：SVD及其应用梁杰存2014310739航博1431.方法求矩阵A奇异值分解一个途径是求解A T A的特征值，但因为舍入误差容易丢掉小奇异值。
因此通常先将矩阵上双对角化，即构造正交阵Q和W,使得Q T AW=B(upper−bidiagnal)。
这一过程可以通过逐次Household变换或逐次Given’s变换完成，还有一种基于待定系数法思想的Lanczos算法。
由于Linpack中SVD算法需要输入上双对角矩阵，本文采用Lanczos 算法实现上双对角化。
1.1.隐式零移位QR法（implicit zero-shift QR）与传统移位QR迭代算法不同，隐式零移位QR算法不进行移位，并且第一步构造右乘Given’s变换矩阵GR(1,2)将上双对角矩阵B的(1,2)位置上的元素12b消零，而不是传统方法中引入一个非零元素21b。
但这一步可能会使原来为零的b12变为非零。
第二步左乘Given’s阵GL(1,2)使得12b为0，但可能会使为零b13变为非零。
与上述步骤类似，将b13变为0后可能会使b23非零。
如下图所示，重复上述步骤最终将恢复为上双对角矩阵，即完成一步隐式零移位QR迭代。
反复迭代，矩阵B将趋近与对角阵阵，对角元即特征值。
图1隐式QR迭代1.2.分而治之(Divide-and-conquer)分而治之算法将上双对角阵B分成有两个互相独立对角块矩阵与另一矩阵之和，即：B=B100B2+b m vvT=Q1Σ1Q1T00Q2Σ1Q2T+b m vv T =Q100Q2(Σ100Σ1+b m uuT)Q1T00Q2T所以矩阵B的特征值与矩阵D+ρu u T的特征值相同，其中D=Σ100Σ1为对角阵，又：det D+ρu u T−λI=det((D−λI)(I+ρD−λ−1u u T))由于D−λI非奇异，则det I+ρD−λ−1u u T=1+ρu T D−λ−1u=1+ρu i2d i−λ=0ni=1在每个d i与d i+1之间分布着一个特征值，可用牛顿法快速找到该特征值。
数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。
求解方法大致可分为直接法和迭代法两大类。
直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。
当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。
如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。
Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。
对正定对称矩阵，采用平方根方法无需选主元。
方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。
实验一一、实验名称：矩阵的LU分解.二、实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b, 并比较这两种方法.三、实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法。
（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组（1）①不选主元高斯消去法求解方程组function [L,U]=gauss1(A,b)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A(:,1:n),-1)+eye(n);U=triu(A);主程序为：Clear；clc；a=ones([84,1])*6;b=ones([83,1]);c=ones([83,1])*8;A=diag(a)+diag(b,1)+diag(c,-1); d=ones([82,1])*15;b=[7;d;14];[L U]=gauss1(A,b);n=length(A);y=ones(n,1);y=L\b;x=ones(n,1);x=U\y解为：x=11.000000000000001.000000000000001.000000000000000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.9999999998835921.000000000232820.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.9999694842845201.000061031430960.9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811 -0.999877927824961 4.99975585192486 -6.99951168894947 16.9990233182979 -30.9980463981918 64.9960918427676 -126.992179871071 256.984344484284 -510.968627937136 1024.93701174855 -2046.87304699420 4096.74218797682 -8190.46875190732 16383.8750076293 -32764.5000305175 65531.0001220701 -131055.000488281 262097.001953124 -524127.007812498 1048001.03125000 -2094975.124999994185857.49999999-8355328.9999999716645128.9999999-33028126.999999965007744.9999998-125821439.000000234866688.999999-402628606.999999536838144.999998可见，这是一个病态方程，从56个跟开始发散。
第一章 绪论(1)1．设0x >,x 的相对误差为δ，求ln x 的误差。
解：近似值*x 的相对误差为*****r e x x e x x δ-=== 而ln x 的误差为()1ln *ln *ln **e x x x e x =-≈ 进而有(ln *)x εδ≈2．设x 的相对误差为2%，求nx 的相对误差。
解：设()n f x x =，则函数的条件数为'()||()p xf x C f x = 又1'()n f x nx-=, 1||n p x nx C n n-⋅∴== 又((*))(*)r p r x n C x εε≈⋅且(*)r e x 为2((*))0.02n r x n ε∴≈3．下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位，试指出它们是几位有效数字：*1 1.1021x =,*20.031x =, *3385.6x =, *456.430x =,*57 1.0.x =⨯ 解：*1 1.1021x =是五位有效数字； *20.031x =是二位有效数字； *3385.6x =是四位有效数字； *456.430x =是五位有效数字； *57 1.0.x =⨯是二位有效数字。
4．利用公式(2.3)求下列各近似值的误差限：(1) ***124x x x ++,(2) ***123x x x ,(3) **24/x x . 其中****1234,,,x x x x 均为第3题所给的数。
解：*41*32*13*34*151()1021()1021()1021()1021()102x x x x x εεεεε-----=⨯=⨯=⨯=⨯=⨯***124***1244333(1)()()()()1111010102221.0510x x x x x x εεεε----++=++=⨯+⨯+⨯=⨯ ***123*********123231132143(2)()()()()1111.10210.031100.031385.610 1.1021385.6102220.215x x x x x x x x x x x x εεεε---=++=⨯⨯⨯+⨯⨯⨯+⨯⨯⨯≈**24****24422*4335(3)(/)()()110.0311056.430102256.43056.43010x x x x x x xεεε---+≈⨯⨯+⨯⨯=⨯=5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为343V R π=则何种函数的条件数为23'4343p R V R R C V R ππ===(*)(*)3(*)r p r r V C R R εεε∴≈=又(*)1r V ε=故度量半径R 时允许的相对误差限为1(*)10.333r R ε=⨯≈6．设028Y =，按递推公式1n n Y Y -=（n=1,2,…）计算到100Y 27.982≈（5位有效数字），试问计算100Y 将有多大误差？解：1n n Y Y -=-10099Y Y ∴=-9998Y Y =9897Y Y =-……10Y Y =-依次代入后，有1000100Y Y =-即1000Y Y =27.982, 100027.982Y Y ∴=-*310001()()(27.982)102Y Y εεε-∴=+=⨯100Y ∴的误差限为31102-⨯。
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
高等数值分析实验一工物研13 成彬彬2004310559一．用CG，Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b作实验中，A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵，S是1043阶的一个稀疏阵A= sprandsym(S,[],0.01,3);检验所生成的矩阵A的特征如下：rank(A-A')=0 %即A=A’，A是对称的；rank(A)=1043 %A满秩cond(A)= 28.5908 %A是一个“好”阵1．CG方法利用CG方法解上面的线性方程组[x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043);结果如下：Iter=35，表示在35步时已经收敛到接近真实xrelres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差绘出A的特征值分布图和收敛曲线：S=svd(A); %绘制特征值分布subplot(211)plot(S);title('Distribution of A''s singular values');;xlabel('n')ylabel('singular values')subplot(212); %绘制收敛曲线semilogy(0:iter,resvec/norm(b),'-o');title('Convergence curve');xlabel('iteration number');ylabel('relative residual');得到如下图象：为了观察CG方法的收敛速度和A的特征值分布的关系，需要改变A的特征值：(1).研究A的最大最小特征值的变化对收敛速度的影响在A的构造过程中，通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数)，可以达到改变A的特征值分布的目的：通过改变rc=0.1，0.0001得到如下两幅图以上三种情况下，由收敛定理2.2.2计算得到的至多叠代次数分别为：48，14和486，由于上实验结果可以看出实际叠代次数都比上限值要小较多。
《数值分析》第二次大作业题目：SOR最优松弛因子选取方法及SOR迭代法的改进内容：1.SOR最优松弛因子选取方法2.SOR迭代法的改进（SSOR迭代法）3.SSOR迭代法的Matlab程序4.举例比较jacobi，Gauss-Seidel，SOR及SSOR 迭代法的收敛速度姓名：合肥工业大学学号：2011班级：信息与计算科学11-1班参考资料：1.《确定SOR最优松弛因子的一个实用算法》李春光等《计算力学学报》2.《数值分析与实验》，薛毅,北京工业大学出版社.3.《数值分析中的迭代法解线性方程组》，马云,科学出版社4.《非线性互补问题的改进超松弛迭代算法》，段班祥等,江西师范大学出版社5.《迭代法解线性方程组的收敛性比较》，郑亚敏,江西科学出版社.一、SOR最优松弛因子选取方法SOR迭代法迭代公式：x(k+1)i=(1-ω)xi+(k) bi-∑aijxjaii⎝j=1ω⎛ i-1(k+1)-j=i+1∑axijn(k)j⎫⎪ (i=1,2,..n.), ⎪⎭1.二分比较法将松弛因子1/2，ω的区间（1,2）进行二分，每个小区间的长度为ω去中间值3/2，按照SOR 迭代法迭代公式，求出跌代次数k,如果k不超过指定的发散常数，则可确定ω的值；否则将（1,2）四等分，每个区间长度为1/4,ω取各分点值，继续迭代，一般地，将1区间（1,2）二分M次，每次二分步长为，ω一次取取各分点值，2M按照SOR迭代法迭代公式，求出跌代次数k,如果k不超过指定的发散常数，则可确定的ω的值，这样总能找到一个不超过指定发散常数ω值。
2.逐步搜索法将1+ω的取值区间（1,2）进行M等分，ω分别取ω的值。
12M-1,1+,...,1+，通过迭代公式依次对同意精度要求求出迭代MMM次数k的值，并从中选出最优松弛因子3.黄金分割法依据黄金分割比的思想，通过计算机主动选取最优松弛因子的近似值，步骤如下a.对（1，2）区间进行第一次0.618的分割，区间边界a1=1,b1=2，在区间(a1,b1)分割出黄金点p1=a1+0.618(b1-a1)，进行SOR迭代法的迭代，求出迭代次数k的值，如果没有超过规定的发散常数，迭代结束，否则做步骤b。
高等数值分析第二次实验作业Tl・构造例子特征值全部在右半平面时，观察基本的Arnold!方法和GMRES方法的数值性态，和相应重新启动算法的收敛性.Answer：(1)构造特征值均在右半平面的矩阵A：根据实Schur分解，构造对角矩阵D由n个块形成，每个对角块具有如下形式，对应一对特征值a - ± ip t£ 2 ⑴5,=lA e丿这样D=diag(Si,S2,S3……SJ矩阵的特征值均分布在右半平面。
生成矩阵A=U T AU,其中U为正交阵，则A矩阵的特征值也均在右半平面。
不妨构造A如下所示:1 -11 12 -22 2n/2 -〃/2n/2 n/2由于选择初值与右端项:xO=zeros(2*N, l)；b=ones(2*N t l)；则生成矩阵A的过程代码如下所示：N=500 %生成八为2N阶A=zeros(2*X)；for a=l：NA(2*a-112*a-1)=a；A(2*a-L.2*a)=-a；A(2*a,2*a-l)=a；A(2*a,2*a) =a；endU = orth(rand(2*N.2*N))；Al = IT 林*U；(2)观察基本的Arnoldi和GMRES方法编写基本的Arnold!函数与基本GMRES函数，具体代码见附录。
Function [x,rm,flag]=Arnoldi(A f b t xO,tol,m)function [x t rm,flag]=GMRES(A t b t xO,tol.m)输入:A为方程组系数矩阵，b为右端项，xO为初值，tol为停机准则，m为人为限制的最大步数。
输出：x为方程的解，nn为残差向虽，flag为解是否收敛的标志。
外程序如下所示： e=le-6； m=700； tic[xA,rmA t flagA]=Arnoldi (A 1, b, xO t e.m)； toe tic[xG,rmG.flagG]=GMRES(Al,b,x0,e,m)； toesubplot(1,2,1)； semi logy(rmA)title(r /\rnoldiE O A 2 u I R ') xlabel (r p u z u 2 Vi E y k*) ylabelC log(| | rk |)') subplot(1,2,2)； semi logy(rmG)title(r GMRESE 6 A 2 u I R ') xlabel (r p u z u 2 Vi E y k*) ylabelC log(| |rk| |)') 得到：waa迭代次数 546526运行时间(s)Arnold"攵敛曲线10?f --------------------- ■------------------- l101 ：\10° : 102101•1-2 o1 1 (-¥--) 60-10200 400 迭代步数k 600 >方法ArnoldiGMRES10°=亘)60200 400 迭代步数k600GMRES 收敛曲线得到两种方法的收敛曲线如上所示，将计算结果整理在下表中:结果讨论：1.从图中可以看出，基本的Arnoldi方法经过546步收敛，基本的GMRES方法经过526步收敛，基本的GMRES收敛速度会略烘于基本的Arnoldi方法。
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=∙∙-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n Tn ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1nTi j i ji j ee α=+⎛⎫ ⎪⎝⎭∑具有,000n j jA -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j AB --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭ 而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此：1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
非线性方程的解法实验1.算法设计与比较问题提出：非线性方程组的求解方法很多，基本的思想是线性化。
不同的方法效果如何，要靠计算的实践来分析、比较。
实验内容：考虑算法(1)牛顿法(2)拟牛顿法分别编写它们的matlab程序。
实验要求：(1)用上述方法，分别计算两个例子。
在达到精度相同的前提下，比较迭代次数、浮点运算次数和CPU时间等。
1.1程序清单为使用flops统计浮点运算次数，使用MATLAB5.3版本%f1.m原函数f1function y=f(x)y(1)=12*x(1)-x(2)^2-4*x(3)-7;y(2)=x(1)^2+10*x(2)-x(3)-8;y(3)=x(2)^3+10*x(3)-8;end%ff1.m原函数f1的雅克比矩阵function y=ff(x)y(1,:)=[12,-2*x(2),-4];y(2,:)=[2*x(1),10,-1];y(3,:)=[0,3*x(2)^2,10];end%f1.m原函数f2function y=f2(x)y(1)=3*x(1)-cos(x(2)*x(3)) -1/2;y(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;y(3)=exp(-x(1)*x(2))+20*x(3)+1/3*(10*pi-3);end%ff2.m原函数f2的雅克比矩阵function y=ff2(x)y(1,:)=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3))];y(2,:)=[2*x(1),-2*81*(x(2)+0.1),cos(x(3))];y(3,:)=[-x(2)*exp(-x(1)*x(2)),-x(1)*exp(-x(1)*x(2)),20]; end%牛顿法(以第一个方程组为例)clear;x0=[0,0,0]';n=10;tol=1e-6;x(:,1)=x0;i=1;u=[1,1,1]';tic;while (norm(u)>tol*norm(x(:,i))&(i<n))A=ff1(x(:,i));b=f1(x(:,i))';u=-A\b;x(:,i+1)=x(:,i)+u;i=i+1;end;x(:,i)iter=i-1t=toc%拟牛顿法(以第一个方程组为例)clear;x0=[0,0,0]';n=10;tol=1e-6;x(:,1)=x0;i=1;p=[1,1,1]';A=ff1(x(:,1));tic;while (norm(p)>tol*norm(x(:,i))&(i<n))x(:,i+1)=x(:,i)-A\f1(x(:,i))';p=x(:,i+1)-x(:,i);q=f1(x(:,i+1))'-f1(x(:,i))';A=A+(q-A*p)*p'/norm(p,2)^2;i=i+1;end;iter=i-1t=tocx(:,i)1.2运行结果1.2.1第一个方程组精确解为*T =(0.886020214719037, 0.796444775323146, 0.749479574122230)x 取最大迭代次数n=5000，相对误差限Tol=1e-6 (1)取()(0)1,1,1x T=牛顿迭代法迭代3次收敛，浮点运算次数为440，每次迭代平均浮点运算次数为147，CPU 耗时t =0(s)拟牛顿法迭代4次收敛，浮点运算次数为1048，每次迭代平均浮点运算次数为262，CPU 耗时t =0(s)(2)取()(0)000x T =，， 牛顿迭代法迭代4次收敛，浮点运算次数为510，每次迭代平均浮点运算次数为128，CPU 耗时t =1.600e-002(s)拟牛顿法迭代6次收敛，浮点运算次数为1493，每次迭代平均浮点运算次数为248，CPU 耗时t =1.50e-002(s)(3)取()(0)50,5050x T=，牛顿迭代法迭代15次收敛，浮点运算次数为2118，每次迭代平均浮点运算次数为141，CPU 耗时t =1.600e-002(s)拟牛顿法迭代338次收敛，浮点运算次数为88454，每次迭代平均浮点运算次数为262，CPU 耗时t =3.100e-002(s)1.2.2第二个方程组精确解为*T =(0.886020214719037, 0.796444775323146, 0.749479574122230)x 取最大迭代次数n=5000，相对误差限Tol=1e-6(1)取()(0)000x T=，， 牛顿迭代法迭代5次收敛，浮点运算次数为776，每次迭代平均浮点运算次数为155.2，CPU 耗时t =0(s)拟牛顿法迭代6次收敛，浮点运算次数为1635，每次迭代平均浮点运算次数为273，CPU 耗时t =0(s)(2)取()(0)888x T=，， 牛顿迭代法迭代9次收敛，浮点运算次数为1519，每次迭代平均浮点运算次数为169，CPU 耗时t =0(s)拟牛顿法迭代21次收敛，浮点运算次数为5924，每次迭代平均浮点运算次数为282，CPU 耗时t =1.600e-002(s)(3)对于离精确解更远的初值（如()(0)101010x T=，，），在计算中会出现奇异或接近奇异的矩阵，计算结果误差很大或计算根本无法进行下去。
1. 解：由于求解Ax b =等价于极小化2Ax b -，相当于极小化泛函：()()1,2x Ax b Ax b ϕ=-- 从任一k x 出发，沿着()x ϕ在k x 点下降最快的方向搜索下一个近似点1k x +，使得()1k x ϕ+在该方向上达到极小值，最速下降方向为：()()T T k k x A Ax b A r ϕ-∇=--=令1,T k k k k k k p A r x x p α+==+，需要寻找合适的k α使得()()11min k k Rx x αϕϕ++∈=，则()()()()()()1,=T T T T k k k k k k k k d x x p p A A x p b p A p r d ϕϕααα+=∇=+-- 令()0d x d ϕα=，可得： ()()()(),,,,k k k k k k k k k Ap r p p Ap Ap Ap Ap α==则()22,0k k d Ap Ap d ϕα=>，因此上式的k α即为所求 于是得到极小化泛函()()1,2x Ax b Ax b ϕ=--的最速下降法： 1) 选取初值0n x R ∈，计算00r b Ax =-2) k=0,1,2,……若k r ε≤，则停止；ε为事先给定的停机场数；否则，k=k+1()()11;,;,;;T k k k k k k k k k k k k k k k p A r p p Ap Ap x x p r r Ap ααα++===+=-2. 解：()()111k k k k AAAx x x r x p A x x α***----=+-=-其中()1p A A α=-，设A 的特征根为120n λλλ≥≥≥> ，则有()11max k i k AAi nx x p x x λ**-≤≤-≤-当120αλ<<时，()1max 1i i np λ≤≤<，因此此方法收敛当()111n αλαλ-=--即12n αλλ=+时，()1max i i n p λ≤≤取最小值11n nλλλλ-+，此时收敛最快3. 解：设x *为方程组Ax b =精确解，k k e x x *=-，则()1,TT Tk k k E e e -=，则原迭代法等价于()110k k I A I E E I βαβ+⎡++-⎤=⎢⎥⎣⎦令()10I A I B I βαβ⎡++-⎤=⎢⎥⎣⎦，则迭代法收敛等价于()1B ρ<，即()1,1i B i nλ<≤≤，令0B I λ-=，即 ()10n I A ββλαλλ⎛⎫+--+-= ⎪⎝⎭ 当0λ≠时，则有10I A ββλαλ⎛⎫+--+= ⎪⎝⎭设120n μμμ≥≥≥> 是A 的特征根，则有()210101112i i i ββλαμλλβαμλβλβαμβ+--+=-+++=<⇔++<+<则有：()()()11112,1211,0,1211,0i i B i ni n ρβαμβββαμββαμ<⇔++<+<≤≤+⇔<-<<≤≤+⇔<-<< 4.5. 证明：反证法。
